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Preface 


I  hope  this  work  will  prove  useful  to  future  sorption  modeling  and  groundwater 
cleanup  efforts.  The  task  is  indeed  ominous,  and  understanding  all  of  its  nuances  is  far 
from  my  grasp.  The  patient  instruction  of  Maj.  Edward  Heyse  and  much  study  yielded 
the  understanding  I  have  gained.  My  indebtedness  in  groundwater  and  contaminant 
transport  topics  is  entirely  to  him.  Dr.  William  P.  Baker  gave  excellent  tutoring  early  on 
in  the  area  of  differential  and  integral  equations.  The  foundation  of  my  limited 
knowledge  of  numerical  analysis  was  laid  by  Maj.  Dave  L.  Coulliette,  but  the  tough  task 
of  raising  me  to  the  level  of  understanding  in  this  area  necessary  for  this  work  was 
expertly  accomplished  by  my  advisor,  Dr.  Kirk  A.  Mathews.  He  posed  the  problem, 
roughly  formulated  the  model,  proposed  the  Volterra  integral  equation  solution  approach, 
and  suggested  the  trapezoidal  algorithm,  not  to  mention  all  the  timely  advice,  help  and 
direction  he  gave  along  the  way.  I  thank  my  wife,  Corina  for  enduring  this  time-rigorous 
process.  I  thank  my  daughter  Kira  for  her  many  smiles,  giggles  and  self  control  in  not 
consuming  this  work,  though  her  appetite  for  paper  products  is  still  quite  strong.  Lastly 
and  most  importantly,  God  has  been  very  gracious  to  me  all  of  my  life.  This  was  surely 
the  case  during  this  work.  Seeing  that  He  “...giveth  to  all  life,  and  breath,  and  all 
things...”  (Acts  17:25),  all  glory  for  any  good  that  comes  of  this  work  is  due  Him,  because 
I  have  only  done  that  which  is  my  duty  to  do  towards  Him  (Luke  17:5-10,  Colossians 
3:22-25).  Without  the  Lord  Jesus  Christ,  we  can  do  nothing  (John  15:5).  I  owe  Him  all. 
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Abstract 


Many  try  modeling  groundwater  contaminant  transport  to  predict  it.  Is  this 
possible  with  rate-limited  processes,  and  under  what  conditions?  On  occasion,  cleanups 
go  slower  than  predicted  (tailing)  and  hazardous  concentrations  reappear  after  cleanup  is 
thought  complete  (rebound).  Rate-limited  transport  is  blamed  by  many.  When  immobile 
water  is  present,  diffusion  from  varied  sizes  and  shapes  of  immobile  regions  can  cause 
varied  rate  limitations  (due  to  varied  diffusion  path  lengths).  Although  known,  most 
modelers  represent  these  varied  rate-limiting  processes  with  a  single  “representative” 
rate-parameter.  This  can  yield  poor  predictions  for  long-term  experiments,  and  the 
parameter  is  generally  time  and  pump-rate  dependent.  This  model  employs  a  distribution 
of  first-order  rate  parameters  to  investigate  the  effects  of  using  a  single  rate-parameter. 
Spatial  effects  are  ignored  by  using  volume-averaged  concentrations  (a  point,  well-mixed 
model)  and  dilutive  pumping  and  rate-limited  transport  are  modeled  to  isolate  rate-limited 
transport  for  study.  A  three-parameter  Gamma  distribution  defines  the  rate  parameter 
continuum.  A  clean  flow  approximation  is  used  extensively,  and  pulsed  pumping  is 
examined  briefly.  An  effective  time  and  pump-rate  dependence  is  seen  in  the  average 
rate.  Long-term  soil  and  contaminant  transport  characteristics  along  with  uptake  history 
or  good  experimentally-derived  initial  contaminant  presence  are  concluded  as  necessary 
for  accurate  predictions. 
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A  POINT  MODEL  OF  AQUIFER  CLEANUP  WITH  A  DISTRIBUTION  OF 


FIRST-ORDER  RATE  PARAMETERS 


I.  Introduction 


Groundwater  Crucial 

To  Mankind.  Our  groundwater  is  a  precious  commodity.  It  provided  drinking 
water  for  53  percent  of  our  nation’s  population  in  1991  (Masters,  1991:104;  Claborn  and 
Rainwater,  1991:1290)  and  is  used  extensively  in  crop  irrigation,  affecting  all  of  our 
diets.  Preventing  and  reducing  the  contamination  of  our  groundwater  and  improving  our 
cleanup  efforts  are  therefore  of  dire  importance  to  all  of  us. 

To  the  USAF.  Department  of  Defense  (DOD)  installations  are  responsible  for  a 
significant  amount  of  accidental  groundwater  contamination,  and  cleaning  up  these  sites 
is,  therefore,  the  responsibility  of  the  DOD.  The  US  Air  Force  alone  will  spend  an 
estimated  $7-10  billion  on  Installation  Restoration  Program  (IRP)  cleanups  from  1992- 
2002,  and  a  “substantial  portion  of  these  costs  will  be  associated  with  groundwater 
contamination  remediation,”  (Adams  and  Viramontes,  1993:1-2;  Vest,  1992).  Seeking  to 
improve  the  effectiveness  and  efficiency  of  cleanup  efforts  is  very  much  an  Air  Force 
problem,  and  one  which  can  be  given  significant  attention  without  apology  to  the  public 
at  large. 


1.1 


Problem  Statement 


At  present,  the  most  common  method  for  cleaning  up  a  contaminated  site  is  to 
sink  extraction  wells  at  various  locations,  pump  out  the  contaminated  groundwater  until 
the  concentration  in  the  ground  is  deemed  safe,  treat  the  extracted  water  to  remove  or 
make  safe  the  contaminant,  and  return  it  to  the  environment  (Adams  and  Viramontes, 
1993).  Models  are  used  to  predict  cleanup  times  and  costs,  but  often  they  fail  at  this  task 
(Travis  and  Doty,  1990:1465).  Pulsed  pumping  is  a  technique  being  increasingly 
investigated  as  an  improvement  to  the  standard  pump  and  treat  technique.  It  involves 
turning  the  pumps  on  and  off  at  calculated  intervals  to  increase  the  average  concentration 
being  extracted  from  wells,  thereby  increasing  treatment  efficiency  (Hartman,  1994).  The 
pump-off  period  will  hereafter  be  called  the  soak  phase.  Improved  predictions  for  pump 
scheduling  could  be  critical  to  the  success  of  this  newer  technique. 

Any  given  failure  of  a  model  to  predict  could  have  many  different  causes.  Two 
behaviors  that  occur  often  and  are  difficult  if  not  impossible  to  predict  are  tailing  and 
rebound  (Olsen  and  Kavanaugh,  1993:44).  Tailing  is  defined  as  a  slower  than  predicted 
decrease  in  the  concentration  of  contaminated  water  streaming  from  the  extraction  well, 
especially  at  later  time  intervals.  This  obviously  involves  longer  than  predicted  cleanup 
times  and  higher  than  predicted  costs.  Rebound  is  defined  as  contaminant  concentration 
rising  significantly  from  the  level  it  was  at  when  pumps  were  turned  off.  On  occasion, 
cleanups  have  been  declared  complete  (with  well  heads  even  removed,  etc.),  only  to 
rebound  years  later  to  unexpectedly  high  (and  hazardous)  contaminant  concentrations. 


1.2 


This  is  frustrating  for  modelers  and  those  tasked  with  cleanup,  but  it  is  more  importantly 
dangerous  for  groundwater  consumers.  The  cause  of  these  two  behaviors  is  not  known 
for  certain,  but  rate-limited  transport  is  frequently  suggested  (Adams  and  Viramontes, 
1993: 1-4).  Present  one-parameter  rate-limited  models  do  predict  some  of  the  long-term 
tailing  and  rebound,  but  not  all  of  it. 

Tailing  and  rebound  both  could  be  generalized  as  inaccurate  predictions  of 
concentration  versus  time.  For  this  work,  inaccurate  is  defined  as  causing  significant 
schedule  deviations,  significant  cost  overruns,  or  the  unpredicted  reappearance  of 
hazardous  concentrations.  The  first  major  question  is:  Are  accurate  predictions  even 
possible?  The  second  is:  If  accurate  predictions  are  possible,  what  conditions  and  tools 
will  facilitate  accurate  predictions  and  are  these  predictions  worth  the  effort? 

There  is  another  baffling  behavior  that  will  be  investigated.  Experiments  to 
determine  “representative”  rate  parameters  have  found  them  to  be  time  and  pump-rate 
dependent  (Brusseau  and  Rao,  1989:56).  This  is  clearly  related  to  the  inadequate 
modeling  of  tailing  and  rebound  described  earlier.  What  is  the  cause  of  these 
dependencies?  Is  there  only  a  single  rate-limiting  process,  but  one  which  is  time  and 
pump-rate  dependent?  Could  the  presence  of  multiple  rate-limiting  processes  cause  this 
behavior? 

Research  Objectives  and  limitations 

The  first  major  objective  of  this  work  is  to  investigate  the  predictability  of  rate- 
limited  contaminant  transport  if  a  range  of  different  rate  limitations  are  in  fact  present. 
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represented  by  a  distribution  of  rate-parameters  (instead  of  the  common,  single 
“representative”  rate-parameter).  To  do  this,  a  model  is  built  to  focus  purely  on  rate- 
limited  transport  and  to  compare  the  behavior  of  a  one-parameter  case  to  that  of  a  case 
assuming  a  continuum  of  rate  parameters,  defined  by  a  continuous  distribution  function. 
Also,  comparisons  will  be  made  between  the  behaviors  observed  for  various  distributions, 
including  a  bimodal.  Assumptions  concerning  every  other  possible  process  will  be  made 
in  order  to  ensure  that  only  rate-limited  transport  and  simple  dilution  cause  the  observed 
behavior. 

The  second  and  final  major  objective  is  to  investigate  the  solutions  obtained  by 
using  a  distribution  of  rate  parameters  in  order  to  determine  whether  the  presence  of 
multiple  rate-limiting  processes  is  a  possible  explanation  for  difficult  to  predict  effects 
such  as  tailing,  rebound  and  the  variability  of  “representative”  rate  parameters  with  time 
and  pump  rate. 

Others  have  used  a  distribution  of  first-order  rate  parameters.  Connaughton  (et  al, 
1993)  modeled  mass  transfer  using  a  multi-site  (multi-rate,  distributed)  model  where  sites 
were  in  parallel,  but  did  not  compare  the  resulting  solution  to  that  of  a  one-parameter 
model.  Heyse  (1994)  modeled  mass  transfer  in  a  continuously-stirred  flow  cell  with  two 
different  distributed  models,  one  like  that  of  Connaughton’ s  and  another  where  the  sites 
were  in  series.  Although  a  one-parameter  model  was  compared  to  the  distributed  models 
in  his  work,  experiments  were  limited  in  duration.  Well  before  both  of  these  researchers, 
Villermaux  (1974)  modeled  mass  transfer  in  a  chromatographic  column  using  what 
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appeared  to  be  a  distributed  model  with  sites  in  parallel,  but  as  with  Connaughton’s  work, 
his  solutions  were  not  compared  to  a  one-parameter  solution. 

In  summary,  the  effects  (especially  long-term  effects)  of  using  multiple  rate 
parameters  instead  of  one  representative  parameter  was  not  clear  in  the  literature  to  this 
author.  The  overall  goal  of  this  work  is  to  investigate  these  effects. 

Overview 

First,  an  overview  of  general  soil  characteristics  important  to  this  work, 
groundwater  contaminant  storage  and  transport  concepts,  and  modeling  will  be  given  in 
the  Background  section.  In  addition  to  a  brief  overview  of  concepts  discussed  and  used 
in  following  sections,  a  theoretical  justification  for  using  multiple  rate  parameters  is 
sought.  In  the  Model  Development  section,  governing  equations  for  two  contaminant 
transport  models  will  be  constructed  based  on  stated  assumptions  for  a  mobile/immobile 
zone  conceptualization  and  for  a  chemical  conceptualization,  using  three  different  rate 
parameter  cases: 

1.  The  unimodal  (distribution  of  rate  parameters)  case. 

2.  The  bimodal  (distribution  of  rate  parameters)  case,  and 

3.  The  one-parameter  case,  (i.e.  the  case  for  most  present  rate-limited  sorption 

models) 
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with  an  overview  and  analysis  of  distribution  function  choices.  In  the  Solution 
Approaches  section,  various  approaches  to  solving  the  three  cases  are  examined, 
including  a  discussion  of  a  clean  flow  approximation,  pump  rate  changes,  a  soaking 
phase,  and  solution  uniqueness.  In  general,  this  analysis  yields  either  a  Volterra  integro- 
differential  governing  equation,  or  a  Volterra  (second  kind)  integral  governing  equation. 

The  Numerical  Methods  section  examines  three  different  classes  of  numerical 
methods  for  both  the  Volterra  integro-differential  type  and  Volterra  (Second  Kind)  type 
governing  equations.  Error  estimation,  pump  rate  changes,  soaking  and  reinitialization 
are  discussed  there,  along  with  a  justification  for  the  method  of  choice.  In  the 
Applications  section,  the  developed  tools  are  used  to  examine  various  solution  behaviors, 
leading  to  the  fulfillment  of  stated  objectives.  Finally,  Conclusions  and 
Recommendations  are  given,  along  with  perceived  areas  of  uncompleted,  related  work. 
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II.  Background 


Soil  Characteristics  and  Groundwater  Contaminant  Transport  Concepts 

Generally,  many  soil  characteristics  can  affect  contaminant  transport,  and  hence, 
cleanup  operations.  Porosity,  n,  is  the  ratio  of  void  volume  to  total  volume.  Saturation, 

S,  is  the  ratio  of  water  volume  to  void  volume.  Volumetric  water  content,  0,  is  the 
product  of  porosity  and  saturation,  or  the  ratio  of  water  volume  to  total  volume. 

Properties  of  soils  are  stated  on  the  basis  of  the  scale  over  which  they  are  possibly  variant 
or  invariant:  microscopic,  macroscopic  and  megascopic.  For  example,  porosity  may 
appear  invariant  on  the  macroscopic  scale  (sometimes  called  the  Darcy  scale),  while 
varying  significantly  on  the  microscopic  and  megascopic  scales.  Finally,  the 
chemical/elemental  constituents,  surface  properties,  and  distribution  of  solids  in  soil  can 
be  important. 

Typically,  aqueous  phase  contaminants  are  easily  cleaned  up  by  the  process  of 
advection.  Advection  is  the  movement  of  contaminant  due  solely  to  the  mean  flow  of 
water  (Domenico  and  Schwartz,  1990:358).  Cleanup  operations  usually  involve  pumping 
groundwater  out  of  contaminated  sites  at  flow  rates  that  allow  uncontaminated  water  from 
surrounding  regions  to  flow  through  the  site  and  evacuate  a  great  majority  of  originally 
present  mobile  water  in  a  relatively  short  time  (Adams  and  Viramontes,  1993).  In 
laboratory  work,  contaminated  water  is  typically  drawn  into  the  experimental  soil  section, 
and  the  concentration  exiting  the  opposite  end  of  the  experimental  section  is  analyzed. 
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Both  in  laboratory  experiments  and  in  field  work,  when  concentration  versus  time  is 
plotted,  breakthrough  curves  (or  BTC’s)  are  said  to  be  generated. 

The  term  breakthrough  is  used  both  in  field  and  lab  work.  This  is  because 
contaminant  concentration  typically  stays  nearly  constant  for  a  period  of  time.  In  the  lab, 
this  time  period  is  the  time  required  for  the  first  contaminated  effluent  to  arrive  at  sensors. 
In  the  field,  this  period  is  that  required  for  cleaner  flow  from  surrounding  uncontaminated 
sites  to  break  through  the  contaminated  region  of  interest,  arriving  at  well-heads/sensors. 
Once  this  period  is  over,  sharp  increases  (in  lab  work)  or  sharp  decreases  (in  field  work) 
in  contaminant  concentration  follow. 

In  addition  to  advection,  dispersion  and  diffusion  within  the  mobile  region  can 
effect  the  shape  of  BTC’s.  Diffusion  involves  the  movement  of  contaminant  opposite 
concentration  gradients,  and  dispersion  involves  the  spreading  of  contaminant  due  to 
local  variations  in  fluid  flow  velocities.  The  effects  of  advection,  dispersion,  and 
diffusion,  at  least  in  the  mobile  region,  are  usually  dominant  in  the  short  term  as 
compared  to  rate-limited  sorption.  This  is  evident  in  the  tailing  which  occurs  when 
models  do  not  account  for  rate-limited  sorption:  advection,  dispersion  and  mobile  region 
diffusion  are  being  modeled  successfully,  but  rate-limited  sorption  does  not  begin  to 
dominate  until  the  tailing  appears  in  the  longer  term. 

As  we  shall  see  in  greater  detail,  a  well-mixed  situation  as  will  be  used,  ignores 
dispersion  and  diffusion  in  the  mobile  region,  and  reduces  advective  processes  to  simple 
dilution.  Because  of  this,  the  short-term  regions  of  BTC’s  generated  by  this  work’s 
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models  are  not  intended  to  be  accurate  and  are  not  important  to  the  work.  It  will  be  the 
longer  term  that  will  be  of  interest,  due  to  rate-limited  transport  domination  there. 

The  contaminants  that  often  pose  real  cleanup  problems  are  those  in  the  solid 
phase  (sorbed  included).  Solids  typically  do  not  move  with  natural  groundwater  flow  or 
during  pumping  operations.  They  continue  to  supply  contaminant  to  the  mobile  region  to 
allow  removal,  but  at  much  slower  rates  due  to  the  following  processes. 

Contaminant  Storage 

Contaminants  are  introduced  into  groundwater  in  a  variety  of  ways  and  on  a 
variety  of  schedules.  Sometimes  exposures  are  short  and  highly  concentrated,  while  other 
times  they  trickle  in  over  long,  continuous  periods  of  low  concentration  exposure. 
Contaminants  can  enter  the  ground  in  various  phases  and  change  phase  as  they  are 
deposited  or  transported.  Although  this  work  is  primarily  focusing  upon  cleanup 
behavior,  adsorption  behavior  during  deposition  for  a  reversible  process  is  clearly  in  view 
also.  It  is  this  process  which  ultimately  determines  valid  initial  contaminant  levels  for 
cleanup  efforts.  For  reversible  processes,  conclusions  about  desorption  behavior  can  be 
equally  applied  to  adsorption  behavior. 

Sorption  and  Immobile  Regions.  In  general,  we  will  investigate  two  possible 
rate-limiting  processes  which  could  cause  some  of  the  inaccuracies  in  present  modeling 
techniques.  Rate-limited  refers  to  processes  which  limit  the  possible  rate  at  which 
contaminant  can  be  removed  from  storage  in  the  ground.  This  work  is  primarily 
concerned  with  contaminant  storage  by  two  basic  processes;  adsorption  and  diffusion 
into  immobile  regions.  Adsorption  is  the  process  by  which  contaminant  bonds  to  solid 
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soil  surfaces,  and  desorption  is  the  breaking  of  these  bonds.  Although  absorption  is 
normally  included  in  the  term,  sorption  in  this  work  will  be  primarily  referring  to 
adsorption.  Adsorption  bonds  fall  into  three  basic  categories,  with  some  subcategorizes 
based  upon  slight  variance  in  bond  type: 


1.  Chemical,  or  chemisorption 

a.  Covalent 

b.  Hydrogen  Bond 

2.  Electrostatic 

a.  Ion-Ion 

b.  Ion-Dipole 

3.  Physical 

a.  Dipole-Dipole/Coulombic 

b.  Keesom  energy 

c.  Dipole-Induced  Dipole/Debye  energy 

d.  Instantaneous  Dipole-Induced  Dipole/London  dispersive 
energy 

(Weber  etal,  1991:501). 

The  strength  of  a  given  adsorption  bond  as  well  as  the  conditions  under  which  it  will  form 
or  be  broken  is  a  function  of  contaminant  type/phase,  sorption  site,  and  temperature 
(Treybal,  1980).  I  suspect  that  pore  water  velocity  can  effect  the  bonding  process  also. 
These  bond  types  are  presented  as  evidence  of  the  fact  that  different  types  and  strengths 
of  adsorption  bonds  are  indeed  possible  for  some  contaminants. 

Immobile  regions  are  defined  as  regions  of  immobile  but  connected  water.  If 
these  immobile  regions  are  indeed  present,  it  is  assumed  that  contaminant  diffuses  into 
these  regions.  Thus,  removal  of  contaminant  from  these  regions  can  be  rate-limited  due 
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to  the  diffusive  process  required.  The  existence  and  importance  of  immobile  regions  to 
many  cases  is  well  established  (Brusseau  and  Rao,  1989). 

Equilibrium  Sorption.  If  water  with  a  given  concentration  of  contaminant  is  kept 
in  stagnant  contact  with  sorptive  surfaces  for  a  long  enough  time,  an  equilibrium 
relationship  is  established.  That  is,  the  amount  of  sorbed  contaminant  eventually  reaches 
a  constant  value,  at  which  time  adsorption  and  desorption  are  in  balance.  The  sorption 
reaction  rate  determines  the  time  required  for  equilibrium  to  be  sufficiently  established. 

If  sorbed  mass  ,  S,  is  plotted  versus  contaminant  concentration,  sorption  isotherms  are  the 
result  (constant  temperature  is  assumed).  Most  isotherms  can  be  expressed  by  the 
following: 

S  =KfC  (2.1) 

(Weber  et  al,  1991:506).  This  form  defines  a  linear  isotherm  when  n=l,  an  approximate 
Langmuir  isotherm  when  n  is  appropriately  less  than  one,  and  a  Freundlich  isotherm 
generally.  If  this  relation  holds  (with  a  constant  n  and  for  both  adsorption  and 
desorption,  the  process  is  said  to  be  singular.  If  isotherms  are  non-singular,  or  hysteretic, 
some  contaminant  bonds  permanently  to  the  soil,  rendering  desorption  impossible 
(usually  due  to  a  chemical  reaction).  Finally,  most  isotherms  are  generated  under 
equilibrium  conditions,  eliminating  rate-limiting  effects.  Isotherms  which  are  linear, 
singular,  reversible,  and  produced  under  equilibrium  conditions  are  said  to  be  ideal.  In 
this  work,  isotherms  are  assumed  rate-limited,  but  otherwise  ideal. 
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Non-equilibrium.  Kinetic,  or  Rate-limited  Sorption.  For  a  given  laboratory  or 
field  situation,  some  conditions  of  flow  for  a  given  contaminant  and  soil  do  not  allow 
time  for  equilibration.  Much  laboratory  and  field  investigation  has  proven  many  cases  to 
be  rate-limited  (Brusseau  and  Rao,  1989:41).  Now  we  will  summarize  how  these 
processes  are  presently  modeled,  and  how  these  models  are  used. 

Model  Construction  and  Usage 

As  was  previously  mentioned,  models  are  typically  built  to  produce  accurate  BTC 
predictions  for  the  purpose  of  cleanup  scheduling,  budgeting,  and  safety  considerations. 
Laboratory  work  also  normally  involves  model  construction  for  improved  understanding 
of  the  behavior  being  investigated.  It  is  not  the  point  of  this  work  to  argue  for  or  against 
the  need  for  more  accurate  modeling.  The  possible  importance  of  it  is  briefly  mentioned 
in  the  introduction  and  in  this  section.  If  it  is  found  necessary  to  build  more  accurate 
models,  this  work  suggests  one  possible  means  of  improvement.  Since  it  is  clear  that 
many  situations  prohibit  the  assumption  of  equilibrium  (Brusseau  and  Rao,  1989:41), 
non-equilibrium  will  be  the  primary  focus  of  this  work.  But  to  better  understand  non¬ 
equilibrium,  a  brief  overview  of  how  some  past  and  present  modelers  have  assumed 
equilibrium  is  given. 

Local  Equilibrium  Assumption  (LEAL  This  assumption  is  made  by  many 
modelers  in  building  their  solute-transport  equation  because  of  its  simplifying  effects 
(Adams  and  Viramontes,  1993).  Specifically,  it  involves  assuming  equilibrium  is  reached 
“locally,”  meaning  in  the  discretized  spatial  regions  of  three  dimensional  transport 
models.  Often,  isotherms  of  the  form  given  in  Equation  2.1  are  used  with  n=\  (linear). 
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Whenever  it  is  valid,  it  greatly  simplifies  modeling  without  any  noticeable  accuracy 
losses. 


In  order  for  the  LEA  to  be  valid,  the  rate  of  the  sorption  process  must  be  fast 
relative  to  the  other  processes  affecting  solute  concentration  (e.g.,  advection, 
hydrodynamic  dispersion)  so  that  equilibrium  may  be  established  between  the 
sorbent  and  the  pore  fluid  (Brusseau  and  Rao,  1989:41). 

Most  if  not  all  contaminants  follow  a  pattern  of  a  short,  initially  fast  period  of 
adsorption  (due  to  fast  sites  which  might  allow  local  equilibrium)  followed  by  an 
extended  period  of  slow  adsorption  (presumably  due  to  slow  sites  for  which  the  LEA  is 
invalid)  (Wu,  1986:725).  If  the  process  is  reversible,  this  pattern  is  duplicated  during  the 
desorption  process.  Fast  and  slow  here  refer  to  the  rate  of  kinetic  sorption  relative  to 
advection,  primarily. 

During  laboratory  experiments,  when  contaminant  is  introduced  into  columns  or 
cells  of  soil  by  sending  it  through  in  the  aqueous  phase,  resulting  BTC’s  are  either 
symmetric  or  asymmetric  (Brusseau  and  Rao,  1989).  Asymmetries  are  almost  always 
caused  by  rate-limited  sorption,  the  rare  exception  to  this  being  when  they  are  due 
exclusively  to  hydrodynamic  dispersion  (Brusseau  and  Rao,  1989).  Rate  limitations  in 
the  laboratory  evidenced  by  asymmetric  BTC’s  are  very  common,  and  with  “tailing”  and 
“rebound”  in  the  field,  attest  to  the  widespread  importance  of  rate-limiting  effects  to 
sorption  modeling. 

The  frequent  reality  of  rate-limited  sorption  means  the  frequent  sacrifice  of 
accuracy  in  sorption  modeling  whenever  the  LEA  is  made.  This  sacrifice  may  or  may  not 
be  called  for,  depending  upon  the  relative  importance  of  sorption  compared  to  other 
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processes  and  upon  the  level  of  accuracy  required.  Valocchi  (1986:1699)  demonstrated 
that  with  inward  radial  flow,  there  is  a  range  outside  of  \yhich  the  LEA  is  valid,  due  to 
slower  fluid  velocities.  If  a  three-dimensional  flow  field  is  used  in  any  given  model, 
knowing  where  the  LEA  could  safely  be  made  would  save  unnecessary  complexity  and 
possibly  computing  time.  If  accurate  sorption  modeling  is  required,  non-equilibrium,  or 
rate-limited  sorption  should  be  the  starting  point  and  the  LEA  made  only  when  and  where 
it  is  carefully  proven  valid.  Many  researchers  are  going  to  great  lengths  to  follow  this 
very  advice,  several  of  which  have  already  been  cited. 

Rate-Limited  Transport.  There  is  not  just  one  accepted  way  to  mathematically  or 
conceptually  deal  with  rate-limited  sorption/transport.  Non-equilibrium  models  can  be 
divided  into  two  basic  categories:  chemical  and  physical  (Brusseau  and  Rao,  1989). 

Chemical.  This  class  of  models  involves  the  hypothesis  that  sorption 
sites  themselves  can  be  rate-limited,  without  assuming  the  presence  of  any  rate-limiting 
diffusion  processes  occurring  in  immobile  regions  (Brusseau  and  Rao,  1989).  Most 
assume  two  classes  of  sorption  site:  one  rate-limited  and  the  other  at  equilibrium  (Nkedi- 
Kizza  et  al,  1984: 1 124).  The  differences  between  sites  are  attributed  to  different 
chemical/bonding  interactions  for  one  predominant  soil  content  or  sorption  surface 
feature  (say  limestone)  versus  another  (say  organic  material).  Here,  the  rate-limited  sites 

are  typically  represented  by  a  single  rate-parameter,  k  ,  as  in  the  following  expression  for 
a  sorbed  mass  balance: 

^=q(i-f)jf,c^(r)-s(r)]  (2.2) 
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(Nkedi-Kizza  et  al,  1984: 1 124,1 129).  Here,  S  is  the  total  sorbed  mass  for  a  given  class  of 
sites  per  total  mass  in  soil,  F  is  the  percentage  of  equilibrium-controlled  sites,  and  is 
the  distribution  coefficient  for  linear  adsorption.  This  rate  parameter  for  chemical  models 
is  not  derived  from  first  principles  using  known  physical  or  chemical  constants  (for  the 
contaminant,  soil  or  solvent)  but  purely  by  curve-fitting  experimental  data  (Brusseau  and 
Rao,  1989).  Tildes  are  used  here  and  in  following  equations  for  dimensional  variables 
because  most  of  this  work  will  be  done  in  non-dimensional  variables.  This  reduces 
presentational  complexity  throughout. 

Physical.  With  this  type  of  model,  mass  transfer  is  hypothesized  to  be 
rate-limited  due  to  the  presence  of  immobile  zones  (Brusseau  and  Rao,  1989).  The 
sorption  sites  in  both  the  mobile  and  immobile  regions  are  assumed  to  be  at  equilibrium, 
but  mass  transfer  from  immobile  to  mobile  zones  is  assumed  to  be  a  diffusive  process 
(Brusseau  and  Rao,  1989).  This  diffusive  process  is  modeled  using  three  basic 
techniques: 

1.  Explicitly  with  Pick’s  law  (second  order), 

2.  Explicitly  using  a  semi-empirical  first-order  mass-transfer  expression  to 
approximate  Fickian  diffusion  (semi-  because  this  first  order  rate  parameter  has 
been  successfully  estimated  from  first  principles  using  the  diffusion  coefficient 
and  an  assumed  geometry),  and 

3.  Implicitly  with  the  use  of  an  effective  or  lumped  dispersion  coefficient  that 
includes  the  effects  of  sink/source  diffusion,  as  well  as  hydrodynamic  dispersion 
and  axial  diffusion 

(Brusseau  and  Rao,  1989:46). 
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Although  surface  diffusion  is  possibly  important,  it  is  not  investigated  in  this  work. 
Instead,  diffusion  through  the  immobile  zone  water  is  focused  upon.  The  various  sizes 
and  shapes  of  immobile  zones  are  typically  modeled  by  a  single  representative  (mean) 
size  and  shape.  This  technique,  at  least  in  the  time  frames  experimentally  tested,  seems  to 
have  stronger  justifieation  than  the  chemical  model  for  a  majority  of  contaminants  and 
soils  (Nkedi-Kizza  et  al,  1982:475). 

One  advantage  of  this  conceptual  approach,  as  previously  mentioned,  is  that  rate 
parameters  can  sometimes  be  reasonably  estimated  using  a  first-principles  approach, 
starting  from  known  contaminant  and  solvent  diffusion  properties  and  assumed 
geometries.  The  real  advantage  here  is  gained  if  laboratory  experiments  (especially  long¬ 
term)  become  unnecessary.  Unfortunately,  the  approaeh  has  proven  inapplicable  to  some 
cases,  and  rate  parameter  estimates  are  not  always  accurate  enough  even  when  this 
conceptual  approach  is  proven  valid  (Brusseau  and  Rao,  1989). 

If  immobile  zones  are  found  to  be  of  importance  for  a  given  eontaminant/soil, 
knowing  their  sizes  and  shapes  would  be  helpful.  Unfortunately,  these  sizes  and  shapes 
are  quite  illusive.  If  their  sizes  and  shapes  were  exactly  known,  the  most  rigorous  way  to 
model  the  rate  limitation  would  be  to  use  a  second-order,  Fickian  diffusion  approach. 

The  real  range  and  distribution  of  mass  transfer  rates  would  be  best  analyzed  with  this 
approaeh  using  measured  aggregate  size  and  shape  distributions.  One  researcher  required 
ten  aggregate  size  classes  (all  spherically  shaped)  for  accurate  modeling  (Cooney  et  al, 
1983). 
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The  first-order  technique  described  above  was  found  to  be  the  most  common 
explicit  technique  used,  but  the  general  argument  for  the  presence  of  many  mass  transfer 
rates  based  purely  upon  the  presence  of  many  different  diffusion  path  lengths/geometries 
holds  for  either  explicit  technique.  At  this  point  it  is  noteworthy  that  the  first-order 
chemical  approach  and  the  first-order  physical  diffusion  approach  are  mathematically 
equivalent,  while  the  physical  meaning  of  variables  is  different  (Nkedi-Kizza  et  al,  1984). 

When  sorption  in  the  immobile  region  is  ignored,  the  first-order  kinetic  model  for 
physical  diffusion  yields  the  following  immobile  zone  mass  balance: 

0™^^^  =  k\c,n{t)- C™(0]-  (2-3) 


Once  a  given  immobile  zone  geometry  is  assumed  (spherical  in  this  case),  the  effective 
diffusion  coefficient  and  geometry  combine  to  yield  the  following  estimated  first-order 
rate-parameter  (termed  the  mass-transfer  coefficient  in  the  literature): 


k'  = 


D  0. 
m  im 


(2.4) 


(Brusseau  and  Rao,  1989:56).  Here,  D  is  the  effective  diffusion  coefficient,  b  is  the  mean 
spherical  aggregate  radius,  and  /  is  the  shape  factor  for  transforming  the  spherical 
diffusion  model  to  the  kinetic  model.  Even  though  this  equation  is  the  product  of  a  first- 
order  kinetic  model  assumption,  it  is  clear  that  both  the  diffusion  coefficient  and  path 
length  associated  with  a  given  geometry  assumption  significantly  affect  mass  transfer 
rate.  Large/thick  aggregates  cause  slow  diffusive  transfer  rates  and  small/thin  ones  cause 
fast  diffusive  transfer  rates. 
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It  is  important  to  note  here  that  the  “representative”  rate  parameter  is  time  and 
pump-rate  dependent  (Rao  et  al,  1980;  Brusseau  and  Rao,  1989:56).  This  work  will 
investigate  this  effect  also.  The  lumped  dispersion  coefficient  technique  is  not  used  in 
this  work,  so  this  completes  our  overview  of  physical  (diffusion)  modeling  techniques. 

Summary 

If  immobile  zones  are  proven  important  or  sorption  sites  with  different  kinetic 
reaction  rates  indeed  exist  for  a  given  case,  it  immediately  follows  that  different  mass 
transfer  rates  do  coexist.  This  work  will  investigate  how  much  variation  in  mass  transfer 
rates  requires  more  than  one  rate  parameter  in  the  long-term  modeling  effort.  It  will  also 
investigate  the  relationship  that  this  variation  has  to  tailing,  rebound,  and  the  time  and 
pump  rate  dependence  of  the  supposed  “representative”  rate  parameter. 
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III.  Model  Development 


Goal/General  Assumptions 

The  goals  of  this  model  need  to  be  very  clear:  investigate  rate-limited  transport, 
assume  the  presence  of  multiple  rate-limiting  processes,  model  these  with  varied  mass 
transfer  rates,  and  minimize  complexity.  Both  chemical  and  physical  diffusion  models 
are  used  because  these  two  models  together  include  the  widest  range  of  contaminants.  If 
this  work  was  only  applicable  to  a  small  minority  of  contaminants,  its  usefulness  would 
be  greatly  diminished. 

With  the  chemical  model,  rate-limited  sites  were  assumed  in  parallel  (rather  than 
in  series  as  some  researchers  have  done:  Heyse,  1994;  Brusseau  and  Rao,  1989:43)  to 
keep  it  similar  to  the  physical  diffusion  model  and  to  minimize  complexity.  The  choice 
of  approach  with  physical  diffusion  is  more  difficult. 

It  is  doubtful  that  second-order  diffusion  effects  are  important  to  the  basic  premise 
of  this  work:  simultaneously  present  rate-limiting  processes  causing  multiple, 
simultaneously  present  mass  transfer  rates.  In  addition,  second  order  Fickian  diffusion  is 
likely  to  inject  significant  complexity  into  the  solution  process.  For  these  two  reasons,  a 
first-order  mass  transfer  rate  for  both  models  was  assumed. 

Because  one-parameter  first-order  models  have  failed  to  model  tailing  and 
rebound  adequately  in  some  cases,  and  because  theory  suggests  the  presence  of  multiple. 
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simultaneously-active  rate-limiting  processes  for  many  cases,  a  continuum  of  possible 
rate  parameters  will  be  used. 


Present  Model-Building  Techniques 

Models  that  consider  spatial  effects  typically  involve  using  a  spatial  discretization 
(an  elementary  volume)  which  assumes  soil  parameters  of  interest  are  somewhat  constant. 
This  elementary  volume  is  normally  shrunk  to  differential  size  for  the  governing 
equations,  becoming  finite  in  size  upon  implementation  of  the  chosen  discretization 
scheme.  All  models  involve  doing  a  mass-balance  on  this  elementary  volume,  setting  the 


time  rate  of  change  of  contaminant  mass  in  the  elementary  volume 


equal  to  the  net 


losses  and  gains  through  the  boundaries  plus  or  minus  any  sources  or  sinks.  This  time 
rate  of  change  of  contaminant  mass  term  is  typically  broken  up  into  all  the  phases  in 
which  a  contaminant  could  be  stored:  possibly  multiple  solid,  liquid,  aqueous  and 
gaseous  forms.  For  this  work,  there  will  be  three  storage  phase  possibilities:  mobile 
aqueous  phase,  immobile  aqueous  phase,  and  adsorbed  solid  phase. 

The  other  side  of  the  mass  balance  (losses/gains  through  the  boundaries  and 
source/sinks)  represents  how  mass  moves  between  elementary  volumes  and  whether  or 
not  mass  could  be  effectively  created  or  destroyed  (an  example  of  this  would  be  if  a 
contaminant  becomes  non-toxic  upon  undergoing  a  certain  reaction,  and  we  are  only 
concerned  about  toxic  forms).  Keep  in  mind  that  this  work  is  not  about  mobile  region 
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contaminant  transport,  but  primarilly  transport  of  contaminant  into  and  out  of  other 
storage  regions.  Sources/sinks  are  also  excluded  in  this  work. 

After  this  mass  balance  is  accomplished,  it  is  usually  restated  in  terms  of 
concentration  units.  An  example  in  the  literature  of  this  sort  of  model  (without  the 
multiple  rate-parameters),  which  considers  mobile  region  contaminant  transport  in  detail 
is: 


dt  dt 


-e  D  ® 

dz 


-e;„v,-^(3.i) 


(Nkedi-Kizza  et  al,  1984: 1 123).  Here,  F’  is  the  proportion  of  sorbed  solid  in  contact  with 
mobile  water,  is  the  apparent  diffusion  coefficient.  Cm  is  a  volume-averaged 
immobile  concentration,  and  Vm  is  the  average  mobile  pore  water  velocity.  This  mass 
balance  assumes  F’,  0  and  p  to  be  constant  with  time  (an  assumption  made  in  this  work 
also),  and  Qm  and  Dm  to  be  independent  of  spatial  position.  The  first  two  terms  on  the 
LHS  represent  the  change  in  contaminant  mass  storage  in  the  mobile  region  with  time 
(allowing  for  aqueous  and  sorbed  phases).  The  second  two  terms  represent  the  change  in 
contaminant  mass  storage  in  the  immobile  regions  with  time  (again,  allowing  both 
aqueous  and  sorbed  phases).  The  terms  on  the  RHS  represent  dispersive/diffusive 
transport  and  advective  transport  respectively.  Keep  in  mind  that  C  here  is  a  function  of 
both  time  and  spatial  position. 
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Overall  Assumptions 


Since  the  main  goal  of  this  model  is  to  single  out  rate-limited  transport  and 
prevent  any  other  process  from  obscuring  the  clear  observation  of  its  effects,  several 
unrealistic  assumptions  are  made.  First,  a  well-mixed  cell  assumption  is  made  to 
eliminate  all  of  the  spatial  dependencies.  This  means  that  once  contaminant  mass 
exits/enters  immobile  regions,  desorbs  from/adsorbs  to  a  solid  surface,  or  enters/exits  the 
volume  of  interest  by  in-flow/out-flow  of  aqueous  contaminant,  it  instantaneously  mixes 
to  cause  an  increase/decrease  in  contaminant  concentration  throughout  the  volume  of 
interest.  This  assumption  allows  the  elementary  volume  to  encapsulate  the  entire  region 
of  interest,  rather  than  just  one  small  spatially-discretized  portion  of  it.  The  mobile 
region  is  therefore  represented  by  one  volume-averaged  concentration  that  is  assumed  to 
be  only  time-dependent.  The  immobile  regions  will  be  represented  as  in  Equation  3.1,  by 
volume  averaged  concentrations. 

Second,  the  boundary  conditions  for  the  region  of  interest,  having  a  mobile  water 
volume  of  Vm,  is  that  mobile  water  is  being  removed  at  a  flow  rate,  Q,  and  water  of  a 
constant  concentration  is  entering  the  region  at  the  same  rate.  This  assumption  yields  the 
following  for  the  RHS  of  Equation  3.1: 

MassFlowIn  —  MassFlowOut  G  ^  ^  ^ 

V  “  w  ®  m  V^O  “  ’  (3*2) 


and  for  clean  water  flowing  in  from  surrounding  regions  (Co=0), 


MassFlowIn  —  MassFlowOut  Q  ^  ^ 
- =  -—^mCm- 

y  y  m  m 

^  ^  m 


(3.3) 
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This  will  represent  contaminant  mass  transport  across  the  boundaries  of  the  region  of 
interest. 

A  First-Order,  Physical  Model 

The  Dimensional.  One-Parameter  Model.  Since  Equation  3.1  is  a  physical 
model,  we  will  start  with  that  type.  The  first-order  rate,  immobile  zone  mass  balance  is: 

ei„^+(i-F')p^  =  r(c„-cJ.  (3.4) 

Notice  that  in  this  form,  the  rate  parameter,  k\  converts  concentration  difference  units  to 
units  of  time  rate  of  change  of  contaminant  mass  stored  in  an  immobile  region  per  total 
volume  of  interest.  The  first  assumption  here  is  that  since  diffusion  through  immobile 
water  is  assumed  to  be  the  main  rate-limiting  process,  there  is  assumed  to  be  no  sorption 
occurring  in  either  mobile  or  immobile  water  (i.e.  ■S'n,=0  and  5im=0).  This  assumption  is 
made  to  ensure  that  equilibrium  type  behavior  is  only  observed  due  to  immobile  zones 
that  could,  in  fact,  be  rate-limiting  given  the  proper  conditions. 

This  does  minimize  overall  contaminant  storage  capacity,  compared  to  sorptive 
possibilities,  but  it  is  the  main  rate-limiting  process  that  we  desire  to  focus  upon.  Note 
that  after  the  completion  of  this  work,  it  was  realized  that  equilibrium  sorption  in  the 
immobile  zones  effectively  lowers  each  rate  parameter,  and  depending  on  the 
contaminant,  could  slow  mass  transfer  significantly.  But,  since  our  non- 
dimensionalization  process  will  be  designed  to  scale  the  model  to  show  behaviors  for  any 
given  average  mass  transfer  rate,  the  initial  ignorance  of  the  effects  of  this  assumption 


3.5 


should  not  effect  any  conclusions  drawn  from  this  model.  It  remains  clear  that 


equilibrium  sorption  in  the  mobile  zone  would  only  serve  to  mask  rate-limiting  effects. 

You  might  recall  from  Equation  2.4  that  fc’ includes  dm,  but  for  our  usage  we 
desire  not  to  include  it.  This  rate  parameter  without  it  will  be  denoted  k.  So,  with  the 
above  stated  assumptions  the  mass  balances  for  immobile  and  mobile  regions  become 
(with  dependencies  shown): 


(3.5) 

Q{t) 

ar  dr  “  v„ 

(3.6) 

This  is  the  model  that  is  to  be  used  for  the  one-parameter  case. 

The  Continuum  of  Rate  Parameters/Immobile  Zones.  The  first  assumption  for  the 
case  of  a  continuum  of  rate  parameters  is  that  the  rate  parameter,  k,  defines  the  class  of 
immobile  zone  being  described  as  well  as  dm.  This  means  that  water  content  is  constant 
for  all  immobile  zones  with  a  given  rate  parameter.  This  is  not  overly  restrictive  in  that 
the  volume  of  immobile  water  is  very  likely  to  have  some  relation  to  the  diffusion 
pathlength,  assuming  the  same  geometry  is  being  compared. 

The  second  assumption  involves  the  definition  of  the  distribution.  If  the  overall 
mass  balance  for  all  immobile  zones  was  to  be  written  as  a  summation,  it  would  be: 


^im)  total 

i=d 


3.6 


Converting  this  to  an  integral  yields: 


Note  that  g  here  has  units  of  0  per  unit  k.  The  bounds  of  integration  here  are  the  bounds 
of  k,  such  that  if  upper  and  lower  limits  of  k  were  to  be  set  differently,  so  would  these 
bounds.  So  to  summarize,  the  overall  mass  balances  for  the  distribution  of  rate 
parameters  (unimodal)  case  become  (with  dependencies  shown): 

e  = e  ™  (*>  [(c„  (<)  -  c, «(*:.'■))] .  (3-9) 


a  f  Q  ot  V  ^ 


(3.10) 


Non-Dimensionalization.  The  following  seven  relations  were  used  to  non- 


dimensionalize  the  model: 


c{t(t))  = 


Cm{t) 


Cim  (^>  f  ) 


t(?)—  kr7. 


m=-, 

kr 


f(m) = -^4— = Ix’ 

^im)  total  ym^r 


(3.11) 
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Note  that  during  this  process,  the  non-dimensional  distribution  is  now  normalized 

(i.e.  f(X)dk  =  1 ).  The  subscript,  r,  denotes  a  constant  reference  value  for  the  given 

variable.  These  reference  values,  when  used  to  non-dimensionalize  this  linear  set  of 
equations,  accomplish  a  very  useful  purpose,  as  the  following  example  demonstrates.  If 
X,  y  and  X  are  assumed  equal  to  one,  solutions  can  be  explored  for  various  values  of  p  and 
V  (minimal  set  of  independent  parameters),  and  then  reference  concentrations  and  rate 
constants  can  be  chosen  to  scale  all  experiments  to  a  given  contaminant  and  set  of  initial 
conditions.  The  appropriate  variable  substitutions  and  chain-rule  derivatives  result  in  the 
following  set  of  governing  equations: 

>  (3-12) 

•  (3-13) 

Dots  above  variables  here  represent  the  first  derivative  with  respect  to  non-dimensional 
time.  These  are  the  two  governing  equations  on  which  this  work  is  based.  The  identical 
non-dimensionalization  for  the  one-parameter  case  yields: 

y(0  =  !^W0-K0]’  (3-14) 

^(0  =  >'(0  ’  (3-15) 

where  |X  here  is  used  for  X  because  in  our  comparisons,  the  mean  of  the  distribution  of  A,’s 
(f(X))  will  be  used  for  this  case.  This  one-parameter  case  will  be  used  not  only  to 
illustrate  the  effects  of  the  one-parameter  assumption,  but  also  as  a  benchmark  problem 
for  testing  numerical  methods.  It  has  an  analytic  solution  and  can  be  solved  by  the  same 
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numerical  techniques  to  be  used  for  the  distribution  of  rate-parameter  cases  which  do  not 
have  analytic  solutions. 


A  Chemical  Model 

Dimensional  Form.  A  chemical  model  will  now  be  built  with  very  similar 
assumptions.  From  the  same  work  as  Equation  3.1,  typical  starting  mobile  water  and 
sorbed  mass  balances  for  this  model  are: 


^  3C„  05.  dSo  ^  ^  ^  ac„ 


(3.16) 


■Sj  =  FK^C^  (equilibrium  sites), 

dt  dt 


(3.17) 

(3.18) 


S2  =  (1  -  F)KjC^  (kinetic  sites  if  equilibrated), 


dS 


dt 


(3.19) 


^  [(1  -  F)K^ (kinetic  site  first  order  rate  law)  (3.20) 


(Nkedi-Kizza  et  al,  1984: 1 124).  Here,  Si  is  sorbed  contaminant  mass  per  total  mass  of 
solids  in  the  volume  of  interest  on  sites  that  instantaneously  equilibrate  and  S2  is  sorbed 
contaminant  mass  per  total  mass  of  solids  in  the  volume  of  interest  on  kinetic  sites.  The 
fraction  of  instantaneously  equilibrating  sites  is  F.  Notice  that  the  first-order  rate 
parameter  here  has  different  units  than  previously,  and  hence  its  symbol  is  different. 

The  same  assumptions  are  made  as  previously  about  the  RHS  of  Equation  3.16:  a 
well-mixed  cell,  volume-averaged  concentrations,  a  total  pumping  rate,  Q,  and  clean 
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water  flowing  in.  With  these  assumptions,  and  upon  substituting  Equation  3.18  into  3.16, 
the  aqueous  phase  mass  balance  can  now  be  written: 

{K+PFKj)^  +  p^  =  -^»„C„.  (3.21) 

dt  dt 

Non-Dimensional  Form.  With  the  following  non-dimensional  transformations 
and  new  definitions  of  v  and  p, 


11 

il-F)K,C,  ’ 

II 

II 

(3.22) 

p(,u 

Vjr{^„  +  pFKj)  ' 

Equations  3.14  and  3.15  are  the  result.  This  demonstrates  for  the  one-parameter  model, 
that  a  certain  set  of  assumptions  can  be  made  to  make  a  chemical  model  mathematically 
equivalent  to  a  physical  first-order  rate  model  in  their  non-dimensional  forms. 
Coefficients  and  variables  have  different  meanings,  but  ultimately  the  actual  behavior 
being  modeled  can  be  equivalently  done  with  either  conceptualization. 

A  Continuum  of  Rate  Parameters.  A  slightly  different  distribution  function 
definition  is  required  for  the  chemical  model.  The  simplest  way  to  incorporate  this 
continuum  of  rate  parameters  is  to  do  so  in  the  non-dimensional  form.  If /(A.)  is  now 
defined  as  the  portion  of  the  total  non-dimensional  sorption  rate  ( )  at  kinetic  sites 

with  non-dimensional  rate-parameter  between  A  and  X+dX  per  unit  A,  Equations  3.12  and 
3. 13  are  the  result.  So,  we  find  at  least  for  this  given  set  of  assumptions  that  these  two 
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approaches,  first-order  chemical  and  physical,  with  an  assumed  continuum  of  mass 
transfer  rates  are  mathematically  equivalent  in  their  non-dimensional  forms.  This 
demonstrates  that  this  model  is  widely  applicable  because  it  accomodates  the  possibility 
of  multiple  mass-transfer  rates  for  either  the  chemical  or  physical  conceptualization. 

Distribution  Choices 

Reality.  Although  theory  and  experimentation  strongly  suggest  the  simultaneous 
presence  of  multiple,  rate-limiting  processes,  a  “representative”  rate  parameter  for  each  of 
these  processes  is  evasive.  If  a  continuum  of  rate-limitations  is,  in  fact,  present,  the 
distribution  is  unknown  for  any  contaminant  for  which  it  is  applicable.  No  one  has 
successfully  measured  the  shape  of  a  given  contaminant’s  rate-parameter  distribution  (as 
far  as  this  author  knows),  although  Heyse  (1995)  is  presently  working  on  this. 

The  diffusion-based  physical  model  suggests  the  usage  of  a  distribution  somehow 
related  to  the  distribution  of  immobile  zones  and  their  associated  diffusion  lengths. 
Lognormal  distributions  are  a  good  experimental  approximation  to  the  sizes  of  aggregates 
in  many  cases  (Gilbert,  1987:152).  It  is  therefore  likely  that  the  lognormal  distribution  is 
somehow  linked  to  the  actual  distribution  of  mass-transfer  rates  for  cases  where  diffusion 
in  and  out  of  aggregates  filled  with  immobile  water  is  a  rate-limiting  process  of 
importance. 

For  the  chemical  model,  multiple  discrete  mass  transfer  rates  makes  better 
theoretical  sense  than  a  continuum.  In  spite  of  this,  conclusions  associated  with  an 
unknown  continuous  distribution  should  be  equally  applicable  to  an  unknown  discrete 
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distribution.  Also,  most  if  not  all  conclusions  associated  with  using  a  single  parameter  to 
represent  multiple  mass  transfer  rates  can  be  based  on  the  behavior  of  models  that  use 
either  a  discrete  or  a  continuous  distribution  of  rates. 

Practicality.  One  general  quality  of  the  distribution  function  to  be  chosen 
determines  its  practicality  for  this  usage:  integrability  in  closed  form  of  various  moments 
of  it.  The  Gamma  distribution  was  found  very  practical,  because  all  desired  integrals 
were  found  analytically,  and  in  closed-form  in  most  every  case. 

Flexibility.  One  additional  criterion  for  this  choice  was  that  the  distribution’s 
shape,  including  its  mean,  standard  deviation,  and  minimum  values,  be  flexible  in 
enabling  desired  variations.  The  three-parameter  Gamma  distribution  has  the  desired 
flexibility  and  closed-form  integrability.  It  also  can  approximate  the  lognormal 
distribution  with  the  appropriate  choice  of  parameters.  Figure  3.1  illustrates  the 
distributions  mentioned. 

The  Gamma  Distribution  (Three-parameter).  This  distribution  is  of  the  form: 


/(>-)= 


1 

a-1 

'-(l-y)' 

pr(a) 

1  p  J 

exp 

[  P  J 

0 

(X>y) 

(^<y) 


(3.23) 


where  -oo  <  y  <  <»,  a  >  0,  P  >  0,  and 


r(a)=  [  (Gamma  function,  normalization  constant).  (3.24) 

■'0 

Sometimes  P  here  is  defined  differently  (as  the  inverse  of  this  parameter:  compare 
Connaughton,  1993:2399  to  Gilbert,  1987: 157).  This  distribution’s  mean  is  |i=aP+y,  and 
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its  standard  deviation  is  (5=^^/a  .  Its  widely  variable  shapes  are  seen  in  Figure  3.2, 
where  the  mean  of  each  distribution  is  (i=l. 

A  Bimodal  Distribution  Allowed.  This  case  only  involves  a  slight  modification  to 
the  variables  found  in  the  governing  equations.  A  new  variable  is  needed  to  specify  the 

percentage  of  the  total  found  in  each  distribution,  di : 

/W  =  <il/lW+(l-dl)/2W'  (3.25) 

Here,/;  and/2  are  Gamma  distributions  combined  in  a  weighted  fashion  (according  to  di, 
with  0<d\<l)  to  obtain  a  normalized  bimodal  distribution.  This  change  retains  the  same 
governing  equations  and  distribution  integrability  while  increasing  distribution  flexibility. 
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IV.  Solution  Approaches 


Note  that  throughout  the  rest  of  this  work,  although  the  chemical  model  yielded 
the  same  non-dimensional  equations,  all  terms  and  discussion  will  favor  the 
mobile/immobile  zone  conceptualization.  The  correlation  between  these  two  types  of 
models  would  be  useful  for  all  discussions  presented  here,  but  time  prevented  the 
inclusion  of  this. 

The  first  step  to  a  solution  is  to  pick  a  solution  approach.  An  approach  that  makes 
sense  and  is  necessary  to  get  to  various  final  forms  is  to  solve  the  immobile  zone 
governing  equation  with  an  integrating  factor.  The  unimodal  case  will  be  handled  first. 

Integrating  Factor  for  Immobile  Zones 

This  technique  is  summarized  with  the  following  four  steps: 

y(?^,r)  =  X.[x(r)-y(X.,r)], 

^  +  y(A.,r)A,  =  x(t)'k, 

at 

^  (4.1) 

-  y{^,0)= 

•'0 

y(A,,r)  =  y(A,,0)e“^^ + (4.2) 

First,  Equation  3.12  is  substituted  into  Equation  3.13,  then  the  above  solution  will  be 
inserted  to  yield  the  following  integro-differential  equation: 
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x(?)= /(^)^  c?A,  .(4.3) 

If  the  order  of  integration  is  switched,  and  we  define  the  following  terms: 

q(t)=  p{t)+\\l, 

(4.4) 

w(r)  =  vj^  f{X)K,y{X,0)e~^‘ dk, 

j(t  -  5)  =  V  J”  / 

our  integro-differential  equation  becomes: 

^(0  =  “^(O-^CO  +  ^(0  +  Jo  ~  •  (4-5) 

With  the  three-parameter  Gamma  distribution  defined  earlier,  and  an  initial  condition  of 
j(X.,0)=>’o  the  integrals  become: 

^  =  Jo°°/(^)^^  =  ap  H-y, 

'^(O  =  '^  Jo~  / e-^‘  dX^yo  , 

.X  f~  .  X  .  .  (4.6) 

;(i;)  =  vJq  f{X)X^e-^^dX, 

v|^|x^-i-a  P^(l  +  2yT)-i-Y^T^^2p 
-  [1+Pt]“-^2 

Recall  that  the  distribution  is  defined  to  be  zero  for  0<  ?i  <  y.  Mathematica  was  used  to 
obtain  these  closed-form  integrals,  using  the  change  of  variable:  z=X-y  (see  Appendix  A). 
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The  bimodal  case  requires  several  different  variable  definitions: 


1^1  = 

^2  “  Jq  f 

m  (0  =  Jo  / 1  Q^)'>^y(^,o)e~^^dk, 

W2{t)  =  vj^  f2{X)Xy{X,0)e-^‘dl,  (4.7) 

Other  than  these  definitions,  the  bimodal  case  is  equivalent  to  the  unimodal  case. 


Integrating  Factor  for  Mobile  Zone 

Generally,  the  integrating  factor  is: 

IntegratingFactor  = 


(4.8) 


but  if  q  is  assumed  constant  with  time,  the  integrating  factor,  with  the  definition  of  W  is: 


IntegratingFactor  =  , 

W{t)=\'  w{s)e^^ds. 


(4.9) 


This  integrating  factor  yields: 


4) = wo)+ «'(')}+ j„vs'w*  - 


\ds. 


(4.10) 


Reversing  the  order  of  integration  and  defining  new  terms  yields: 
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(4.11) 


x'  =  t-s, 

X  =t-s', 

jc(f)  =  e-9'{x(0)+ W(r)}+ (4.12) 

This  form  could  easily  be  altered  to  match  the  general  form  of  a  linear  Volterra  integral 
equation  of  the  second  kind,  and  it  also  includes  a  difference  kernel.  The  importance  of 
this  is  discussed  in  the  solution  uniqueness  section.  Allowing  p  to  change  with  time 
makes  calculating  both  proper  forms  of  W  and  J  a  numerical  endeavor  instead  of  closed- 
form  integrals. 

Analytic  Solution  for  One-Parameter  Model 

The  one-parameter  case,  Equations  3.14  and  3.15,  has  an  analytic  solution  given 
by  (see  Appendix  B  for  the  details): 

z  =  ^(\i-pf  +  pv(2p-H2p-i-|xv), 

Cl  =  —  [xo (-(t  +  P  +  Itv  -I- z)  -  2yQ  pv ], 

C2  =  ^  [xo  (l^  -  P  -  I^V  -t-  z)  +  2yo  pv  ], 

C3  =  —[yo{\^-p-\^^+z)-2xo[^\  (4.13) 

C4  =  +  P -H  pv  -I-  z)  +  2xop]> 


ji:(t)=  ciexp 

-(p  -1-  p  +  pv  -1-  z)t 

+  C2exp 

-(p-l-p  +  pv -z)t 

2 

2 

y(t)=c3exp 

-(p  -1-  p  -1-  pv  +  z)t 

+  C4exp 

-(p-t-p-i-pv-z)t' 

2 

2 
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where  |i  is  defined  in  Equation  4.4.  It  must  be  remembered  at  this  time  that  if  |i  is  a 
function  of  time,  (as  some  researchers  claim:  Brusseau  and  Rao,  1989:56;  Rao  et 
al, 1980a;  Cussler,  1984),  this  approach  would  have  to  be  modified  by  changing  the 
integrating  factor  to: 

IntegratingFactor  =  gl  .  (4. 14) 

Initial  conditions  must  be  reset  anytime  p  is  to  be  changed  to  a  new  constant  value.  This 
solution  was  computed  in  the  main  code  (Appendix  G),  and  is  algebraically  arranged  this 
way  to  minimize  calculations  within  an  iterative  loop.  The  constants  and  even  the 
coefficient  of  t  in  each  exponential  could  be  calculated  outside  the  loop,  and  would  only 
have  to  be  recalculated  at  pump  changing  times. 


Clean  Flow  Approximation 

Another  solution  approach  chosen  assumed  the  pump  rate  was  high  enough  to 
ensure  that  x«y  for  at  least  some  amount  of  time.  This  physically  means  that  clean  water 
is  assumed  to  be  passing  all  immobile  zones  to  determine  mass-transfer  rates  (the  best 
possible  mass  transfer  rate  from  immobile  zones).  Under  these  assumptions,  the 
governing  equations  and  corresponding  solutions  are  (assuming  constant  pump  rate): 

i:(t)  =  -px(t)-vJ“/(>.){-^y(X,r)}dX,  (4.15) 


-pt 


x(0)+  f  wis)e^^ds 

•^0 


(4.16) 
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Here,  w  is  the  function  defined  earlier  in  Equation  4.4.  Equation  3.23  defines  a  to  be 
positive,  and  if  it  is  an  integer,  this  integral  can  be  evaluated  by  recursive  integration  by 
parts,  using  the  relation: 


-1 

m-1  a:'”"* 
until  the  final  integral  is  of  the  form: 


(4.17) 


El(z)=r^s.  (4.18) 

J-Z  s 

Mathematica  was  unable  to  accomplish  this  process,  so  it  was  completed  by  hand  and  a 
Mathematica  module  was  written  to  accomplish  it  (using  the  built-in  ExpIntegralEi[2], 
see  Appendix  C  for  the  details  on  this  module  and  the  function  defining  the  clean  flow 
approximation).  Non-integer  values  for  a  are  assumed  to  yield  a  final  integral  that 
requires  numerical  approximation  (just  as  El  does,  though  a  built-in  Mathematica 
function  would  be  unavailable). 


Pump  Rate  Changes/Soaking 

Here,  problems  associated  with  pump  rate  changes  and  soaking,  specific  to 
solution  approaches  will  be  discussed.  Problems  specific  to  numerical  methods  will  be 
discussed  in  the  Numerical  Methods  section. 

The  General  Integrating  Factor  for  the  Mobile  Region.  The  soaking  phase  is  a 
definite  requirement  in  this  work  to  show  rebound  differences  between  different  models, 
and  it  had  to  immediately  follow  a  given  period  of  pumping  (requiring  a  mid-experiment 
change  in  pump  rate).  It  was  also  a  requirement  to  demonstrate  pulsed  pumping 


4.6 


concepts.  As  was  mentioned  previously  in  this  section,  any  solution  utilizing  an 
integrating  factor  for  the  mobile  region,  (solution  approaches  yielding  the  Volterra 
integral  equation  and  limiting  pump-rate  solution)  assumed  a  constant  pump  rate. 

Clearly,  a  soaking  phase  is  incompatible  with  the  clean  flow  approximation.  (A  pump 
rate  of  zero  is  clearly  below  the  threshold  value  for  a  valid  clean  flow  approximation!) 

It  was  briefly  pointed  out  that  a  different  integrating  factor  than  was  used  would 
be  required  if  pump  rate  was  a  function  of  time.  If  this  general  form  for  the  integrating 
factor  (Equation  4.8)  was  used  to  accommodate  pump-rate  changes,  closed-form 
integrations  would  be  impossible  for  these  solution  approaches,  leaving  only  numerical 
options  for  these  integrations.  Avoiding  numerical  approximations  enhances  accuracy 
and  usually  reduces  computation  cost.  Therefore,  close-form  integrations  are  sought 
when  possible. 

The  One-Parameter  Analytic  Solution.  This  solution  also  assumed  a  constant 
pump  rate,  so  pump  rate  changes  had  to  be  dealt  with  as  initial  conditions.  The  solution 
was  valid  until  the  pump  change  time,  then  the  pumps  were  assumed  to  change 
instantaneously.  The  x  and  y  values  at  pump  change  time  were  then  used  as  initial 
conditions  to  a  new  solution,  with  a  new  “constant”  pump  rate.  This  process  of  starting 
with  new  initial  conditions  is  termed  reinitialization.  There  are  problems  associated  with 
this  process  for  numerical  methods,  to  be  discussed  in  that  section. 


4.7 


Solution  Uniqueness 

To  investigate  solution  uniqueness,  the  solution  approach  involving  an  integrating 
factor  for  the  mobile  region  will  be  examined.  This  approach  led  to  equations  4. 12  and 
4.16,  which  are  easily  manipulated  into  the  general  form  of  a  linear  Volterra  integral 
equation  (second  kind).  In  Yosida’s  notation  (1991:145),  the  general  form  is: 

f{s)  =  (p(:f)-  ,  (4.19) 

where  A,  is  a  constant  and  not  the  rate  parameter, /is  a  known  function,  and  tp  is  the 
function  to  be  found.  Yosida  and  others  have  proven  that  this  general  form  has 
“essentially  only  one  solution”  and  that  solution,  when  examined  based  upon  a  series  of 
iterated  kernels  is  “almost  uniformly  convergent.”  (Tricomi,  1985:10-15;  Yosida, 
1991:145-147;  Linz,  1985:29-35)  The  fact  that  ours  also  includes  a  difference  kernel  also 
contributes  to  uniqueness.  Essentially,  this  form  gives  a  great  deal  of  assurance  to  the 
existence  of  a  unique,  bounded,  non-zero  solution.  This  assurance  is  applicable  to  any 
forms  preceding  it  (i.e.  other  approaches  to  the  same  system  of  differential  equations). 
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V.  Numerical  Methods 


Integro-Differential  Equation 

The  first  numerical  method  attempted  was  applied  to  the  Volterra  integro- 
differential  equation,  Equation  4.5.  A  short  experiment  was  conducted  with  forward  and 
backward  Euler  methods  in  order  to  facilitate  the  construction  of  a  trapezoid  integration 
method.  The  program  built  using  this  method  became  the  workhorse  for  much  of  the 
analysis. 

Trapezoid  Method.  This  implicit  method  is  best  summarized  by  the  following 
sequence,  which  transforms  the  integro-differential  equation  into  simply  an  integral 

equation  (not  using  the  integrating  factor): 

q  =  p+V|i, 
z{t)= 

r  ,  (5.1) 

H\t,  x(t),  z(t)j  =  -qx{t)  +  ) + z{t), 

x(t)  =  /f{t,jr(t),z(t)}. 

First,  h  will  be  defined  as  a  non-dimensional  time  step,  and  Z  and  X  will  be  defined  as 
approximations  of  z  and  x.  After  integrating  over  one  time  step,  tn-i  to  4,  we  obtain: 

+  (5.2) 

t  rt-1 

With  initial  conditions  of  Zo=0  and  Xo=xo,  both  z  and  the  integral  in  Equation  5.2  may  be 
replaced  by  a  trapezoid  approximation: 
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(5.3) 


h  h 

2n  ~  ~Z  ^(^o)^(^O  >  ^n)  , ?/j)  +  “  ^{tn)^{.tn  •>  tn)^ 

^  i=l  ^ 

h 

Xn  ~  Xn-l  ^Xn-l  >Z;i-l)+  >  Xn  »Zn)} 

n  =  1,2,..., 

(Linz,  1985: 178).  This  is  easily  seen  to  be  two  equations  in  two  unknowns,  (Zn  and  Xn) 
with  an  implicit,  algebraic  solution.  Although  this  technique  has  third-order  local 
truncation  error  on  one  step,  Z’s  composite  format  yields  a  second-order  global  truncation 
error  (Isaacson  and  Keller,  1994:316;  Burden  and  Faires,  1985:164).  This  was  proven  for 
the  one-parameter  case:  error  was  proportional  to  h  (see  Appendix  D  for  details  of  this). 
It  was  demonstrated  for  the  unimodal  case  using  uncertainty  estimation  in  the  same 
appendix. 

Multistep  Methods.  Multistep  methods  follow  the  previous  general  procedure, 
but  use  higher  order,  interpolatory  integration  schemes  to  approximate  both  integrals 
(Linz,  1985: 178).  Multiple  starting  values  are  required  for  these  methods,  and  numerical 
instability  is  possible  (Linz,  1985: 178).  Though  this  technique  was  attempted  using 
starting  values  for  both  the  one-parameter  and  one-distribution  cases  from  the  analytic 
one-parameter  case  solution,  no  working  code  was  finished. 

Block-Bv-Block  Methods.  This  technique  involves  initially  the  same  procedure 
as  multistep  methods,  including  the  usage  of  higher  order  integration  schemes.  However, 
it  uses  intermediate  points  (between  nh  and  (n+l)/i)  and  some  form  of  interpolation  to 
generate  multiple  equations  and  multiple  unknowns.  These  unknowns  are  typically  found 
in  a  block  of  several  by  solving  the  equations  simultaneously,  hence  the  name. 
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An  example  of  this  method  would  be  to  use  Simpson’s  rule  and  quadratic 
interpolation  to  generate  two  equations  in  two  unknowns.  Linz  (1985:186)  claimed  that 
this  was  easily  proven  to  be  a  fourth-order  technique.  Fourth  order  accuracy  was  obtained 
with  an  easy  test  case  having  only  one  integration  involved,  but  the  interpolation  scheme 
inserted  no  error  due  to  the  solution  being  a  constant:  y(x)=l  (Linz,  1985: 1 18-120). 
Although  the  integration  scheme  is  locally  fifth-order,  and  globally  fourth  in  its 
composite  format,  the  quadratic  interpolation  is  only  third-order  (see  Appendix  E  for  this 
proof).  I  claim  that  the  technique  is  globally  second-order  due  to  the  usage  of  the 
quadratic  interpolation  formula  twice  on  each  step.  This  was  demonstrated  with  a  simple 
test  case,  and  not  the  actual  model  (see  Appendix  F).  A  higher  order  method  could  no 
doubt  be  designed,  but  the  accuracy  was  not  necessary.  Because  the  efforts  involving 
quadratic  interpolation  only  gave  promise  of  second-order  accuracy,  the  same  order 
available  already  in  the  previously  coded  trapezoid  scheme,  these  efforts  were  abandoned 
before  ever  achieving  higher  order  accuracy  with  a  more  involved  scheme. 

Pump  Rate  Changes/Soaking/Reinitialization.  In  solving  the  integro-differential 
equation,  the  only  reinitialization  required  was  with  the  analytic  solution  of  the  one- 
parameter  case.  The  accuracy  of  this  process  was  quite  important  as  it  was  being 
compared  to  the  same  case  being  solved  numerically  to  closely  monitor  error  growth.  For 
output  purposes  only,  y(k,t)  was  found  at  discrete  )i‘s.  This  gave  important  information 
about  the  soak  phase  redistribution  of  contaminant. 
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Volterra  Integral  Equation  fSecond  Kind) 


This  approach  has  promise  of  better  accuracy  for  a  given  step-size  because  it 
appears  to  involve  the  numeric  approximation  of  only  one  integral  instead  of  two. 
Although  the  integrals  associated  with  it  (when  the  Gamma  distribution  is  used)  are  a  bit 
more  cumbersome,  computationally,  and  involve  approximation  of  exponential  integrals, 
the  method,  without  other  complications,  could  still  be  expected  to  yield  better  accuracy 
per  step.  This  method  is  made  much  more  tractable  by  choosing  a  (Gamma  distribution 
shape  factor)  to  be  an  integer,  but  this  is  of  minimal  consequence;  non-integer  values 
were  not  really  necessary  for  this  work. 

The  main  disadvantages  of  this  method,  associated  mostly  with  pump  changes, 
have  been  alluded  to  in  the  previous  section.  It  was  previously  pointed  out  that  any 
method  using  the  general  form  of  integrating  factor  for  the  mobile  zone  (Equation  4.8) 
required  the  numerical  evaluation  of  all  integrations  otherwise  feasible  in  closed-form. 
This  is  only  one  possible  solution  to  pump  changes  for  this  method.  Another  way  to 
accommodate  pump  changes  is  to  reinitialize  each  time  the  pump  is  to  be  changed,  as  was 
discussed  for  the  one-parameter  case.  Unfortunately,  this  too  forces  numerical  evaluation 
of  otherwise  closed-form  integrations,  unless  an  integrable  (closed-form)  interpolating 
function  could  be  found  for  y(A,,0)  at  each  reinitialization.  I  claim  that  these  problems 
associated  with  pump  changes  would  minimize  accuracy  benefits  of  this  method.  In  fact, 
total  error  could  be  worse,  because  more  actual  integral  approximations  would  be 
required  than  for  the  integro-differential  equation,  after  a  pump  change. 
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Because  the  general  approach  to  integral  equations  is  nearly  identical  to  integro- 
differential  equations,  the  applicable  numerical  methods  fall  into  the  same  categories: 
trapezoid,  multistep,  and  block-by-block  methods.  Although  none  of  these  methods  were 
ultimately  applied  to  the  Volterra  equation,  some  of  the  complexity  was  faced  in 
obtaining  the  clean  flow  approximation.  The  same  integral  was  accomplished  there. 

Justification  for  Choice  of  Method 

As  was  previously  mentioned,  pump  rate  changes  had  to  be  dealt  with  carefully. 
Any  method  attempted  was  faced  with  a  jump  discontinuity.  The  trapezoid  method  as 
applied  to  the  integro-differential  equation  maintained  stability  as  long  as  the  pump  rate 
was  constant  on  a  given  step.  This  involved  the  assumption  that  the  rate  change  was 
instantaneous  at  the  end  of  a  step  such  that  the  old  value  for  the  next  step  was  just  the 
previous  value,  but  the  pump  rates  (old  and  new  values  for  the  new  step)  were  now  at  the 
new  value.  This  was  clearly  the  easiest  method  to  deal  with  pump  rate  changes. 

Because  soaking  after  a  pumping  phase  involved  a  pump  rate  change  (from  non¬ 
zero  to  zero),  soaking  was  clearly  easier  with  the  aforementioned  trapezoid  method.  The 
clean  flow  approximation  was  clearly  invalid  for  a  soaking  phase,  because  a  pump  rate  of 
less  than  10  was  found  to  invalidate  it,  and  pump  rate  was  assumed  constant  for  it  (as 
previously  mentioned).  It  was  therefore  useless  for  pulsed  pumping  also.  Although 
dealing  with  the  Volterra  held  initial  promise  of  better  accuracy  per  step,  the  complexity 
was  a  hindrance.  Ultimately,  due  to  time  constraints  and  previously  stated  reasons,  the 
trapezoid  method  was  applied  to  the  integro-differential  equation  as  the  overall  method  of 
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choice.  However,  much  was  learned  and  demonstrated  with  the  clean  flow 
approximation.  If  accuracy  ever  became  more  critically  important,  a  higher  order  method 
would  probably  be  applied  to  the  integro-differential  equation.  In  that  case,  solving  the 
complications  associated  with  pump  rate  changes  would  be  challenging. 

Error  Estimation 

There  are  two  approaches  typically  used  to  estimate  or  bound  error  values,  and  a 
third  was  devised  for  this  work.  The  first  is  the  rigorous  theoretical  approach  which 
delves  into  the  heart  of  numerical  method  design  and  carries  every  possible  source  of 
error  to  the  final  result.  Such  an  approach  is  intractable  for  our  system.  The  second  is 
simply  to  compute  the  same  sequence  of  values  using  a  sufficiently  smaller  step  size,  and 
observe  the  change  in  the  solution.  This  makes  sense  because  it  is  a  very  safe  assumption 
that  the  error  for  the  first  computation  is  going  to  be  orders  of  magnitude  larger  than  the 
second.  This  is  what  defines  a  sufficiently  smaller  step  size.  This  equates  to  the 
assumption  that  the  difference  between  the  two  computations  very  well  approximates  the 
true  error  of  the  first  computation.  This  method  was  applied  in  this  work. 

It  was  also  noted  in  the  model  development  section  that  the  one-parameter  case, 
although  it  had  an  analytic  solution,  could  be  solved  numerically  by  the  same  methods 
applied  to  the  one-distribution  and  bimodal  cases.  As  a  third  approach,  I  hypothesize  that 
the  exact  error  in  the  one-parameter  case,  which  can  be  easily  calculated,  is  a  good 
benchmark  test  of  error  growth  for  two  basic  reasons.  First,  round-off  errors  and  integral 
approximation  errors  are  very  nearly  identical  in  theory,  since  all  the  cases  when  solved 
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numerically  involve  the  same  numerical  method  with  only  minor  differences  in  their 
defining  functions.  Second,  at  least  for  the  earlier  time  frames,  all  cases  yield  very 
similar  numerical  results. 

Even  though  this  technique  involved  nearly  doubling  the  numerical  computations, 
it  fit  very  efficiently  into  the  algorithm  and  proved  a  very  useful  addition  to  the  program. 
Throughout  the  applications  section,  error  estimates  will  be  referenced  as  originating  with 
this  technique  or  the  second  mentioned  above.  The  second  method  mentioned  (using  a 
significantly  reduced  step-size  to  estimate  uncertainty)  is  probably  better  at  absolute  error 
estimation  than  the  third  (assuming  error  in  one-parameter  case  approximates  error  in 
distributed  cases),  so  it  will  be  used  as  much  as  possible. 
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VI.  Applications 


Uncertainty  of  Numerical  Solutions  Displayed 

All  numerical  solutions  displayed  haye  a  standard  criterion  applied  to  them 
concerning  uncertainty.  Unless  otherwise  noted,  step  size  was  chosen  to  render  estimated 
uncertainty  negligible  (fiye  orders  of  magnitude  smaller  than  the  solution  itself  in  most 
cases).  This  estimate  was  obtained  by  the  numerical  solution  of  the  one-parameter  case, 
for  which  exact  error  was  calculated.  Typically,  ^=0.01  was  sufficient  to  accomplish  this 
uncertainty  order  of  magnitude  goal,  so  deyiations  from  this  step  size  will  be  noted.  On 
seyeral  occasions,  this  estimate  of  uncertainty  was  tested  with  the  smaller  step  size  run  to 
estimate  the  decimal  place  of  change  (hence  the  estimated  decimal  place  of  uncertainty, 
as  preyiously  discussed  also). 

Usage  of  XQ,  Distributions,  and  Distribution  Mean,  p, 

In  this  work,  xq,  yo,  and  |i  were  kept  at  unity  for  solution  purposes.  Reference 

yalues  may  be  chosen  in  order  to  scale  the  solutions  to  a  wide  range  of  dimensional,  field 
and  laboratory  scenarios.  Also,  x  and  y  represent  fractions  of  the  reference  concentration. 
The  solutions  inyestigated  assume  equilibrium  has  been  established  between  x  and  y 
before  the  simulations  begin  (they  are  equal  at  unity).  The  model  could  haye  easily  been 
used  to  inyestigate  soaking  adsorption  (or  storage  in  immobile  regions)  by  starting  y  at 
zero,  or  to  inyestigate  time  required  during  soaking  for  equilibrium  to  be  established  by 
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starting  x  at  zero  and  y  at  unity.  No  contaminated  fluid  was  allowed  to  enter  the  region  of 
interest  from  surrounding  regions,  but  retaining  this  possibility  in  this  model  is  quite 
simple. 


Since  X  was  defined  during  non-dimensionalization  by  Equation  3.1 1  as  a  ratio  of 
the  actual  rate  parameter  to  a  reference  value,  a  unity  value  of  |X  was  used  throughout  this 
work.  In  summary,  for  all  solutions  displayed,  the  initial  conditions  were  ;c=l  and  y=l 
and  a  distribution  mean  of  |i=l  was  used. 

Since  the  three-parameter  Gamma  distribution,  with  parameters  of: 


a  =  3, 

P  =  l/3, 

Y  =0, 

|i  =  ap-i-Y  =1, 


(6.1) 


produced  a  distribution  quite  similar  to  the  lognormal  distribution,  these  were  used, 
unless  otherwise  noted. 


Variations  in  v  and  p 

The  first  observation  to  be  made  about  these  two  adjustable  parameters  is  that  they 
have  similar  effects  upon  the  solution  after  rate-limited  transport  begins  to  dominate 
(after  the  break  in  logarithmic  slope):  that  is,  increased  pump  rates  do  something  similar 
to  decreased  v ’s  as  illustrated  in  Figure  6.1.  Ultimately,  neither  seems  to  affect  the 
logarithmic  slope  after  rate-limited  transport  begins  to  dominate,  but  both  affect  the 
concentration  at  which  this  rate-limited  transport  domination  begins  to  occur.  Pump  rate 
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;,=  10,  V=0.1 
p=40,  v=0.1 
p=\0,  V=0.025 


Figure  6.1.  Effects  of  Changes  in  v  and  p  (One  Parameter  Solution) 


increases  cause  rate-limited  transport  domination  earlier  and  at  lower  concentrations,  and 
an  obviously  increased  logarithmic  slope  during  advection  domination.  This  last 
observation  is  true  even  when  there  are  no  immobile  regions  because  advection  is  really 
only  causing  dilution  due  to  perfect  mixing,  as  mentioned  earlier.  Decreases  in  v  cause  x 
to  follow  the  same  dilution  slope  during  advection  domination  (due  exclusively  to  pump 
rate),  but  follow  it  longer  to  a  lower  concentration. 

In  summary,  advection  (dilution)  is  dominated  by  rate-limited  transport  at  the 
level-off  point,  the  time  and  concentration  for  this  point  being  affected  both  by  p  and  v. 
The  logarithmic  slope  after  rate-limited  transport  domination  begins  is  constant  for  this 
one-parameter  case.  But,  as  we  will  see  later  with  more  clarity,  this  is  only  true  for  the 
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one-parameter  case  and  possibly  a  certain  class  of  distributions  to  be  investigated  (with 
non-zero  y).  In  any  case,  at  least  the  starting  logarithmic  slope  seems  to  be  primarily  a 
function  of  the  distribution  case,  somewhat  independent  of  p  and  v. 

After  examining  both  field  cind  laboratory  data,  and  conferring  with  an  expert  in 
the  field  of  groundwater  hydrogeology,  (Heyse,  1995),  I  observed  that  most  BTC’s  have 
starting  regions  where  advection  dominates.  Therefore,  very  few  solutions  were 
investigated  for  pump  rates  and  v  ‘s  for  which  there  is  no  advection  domination  period. 
The  trapezoidal  code  in  Appendix  G  calculates  this  period  for  the  analytic  solution  so  that 
the  output  scheme  can  sufficiently  follow  the  curvature  at  this  point.  The  point  was 
simply  calculated  by  setting  both  exponential  terms  in  Equation  4.13  equal  to  each  other, 
and  solving  for  the  time  which  makes  the  equality  valid.  A  negative  value  output 
signifies  that  advection  never  dominates. 

These  observations  determined  the  minimum  pump  rate  to  be  investigated.  As 
mentioned  previously,  a  period  of  advection  domination  is  common,  so  the  pump  rate  at 


Table  6.1.  Non-Dimensional  Pump  Rates,/?,  and  v  Combinations 
for  No  Advection  Domination 


Water  Content  Ratio,  v 

Non-Dimensional  pump  rate,  p 

1.0 

2.0 

0.5 

1.5 

0.1 

1.1 

0.01 

1.01 

0.001 

1.001 
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which  this  period  became  minimal  was  used  as  the  minimum  pump  rate.  Table  6. 1  gives 
selected  values  of  p  and  v  for  which  there  is  no  advection  domination.  It  appears  from 
this  table  that  an  approximate  rule  of  thumb  for  no  advection  domination  (using  the  given 
initial  conditions  and  distribution)  is:  p  <  1  +  v  . 

If  it  is  known  that  advection  does  not  initially  dominate  for  a  given  scenario, 
smaller  pump  rates  should  be  investigated.  Heyse  (1995)  estimated  that  a  typical  value 
for  V  is  0.1.  Ball  and  Roberts  (1991:1241)  have  estimated  from  0.004  to  0.2  for  a  similar 
constant.  For  this  work,  unless  otherwise  noted,  a  v  of  0. 1  will  be  used. 

Because  the  clean  flow  approximation  was  easy  to  determine,  its  use  will  be 
maximized.  This  equates  to  using  it  at  all  pump  rates  greater  than  that  for  which  it  is 
deemed  adequate.  Because  the  above  analysis  suggests  that  v  affects  the  value  for  which 
the  clean  flow  approximation  is  adequate,  using  a  value  of  v  different  from  0.1  would 
require  a  repeat  of  the  adequacy  analysis.  For  the  numerical  solution,  a  lower  bound  for  p 
was  set  at  1.0  and  an  upper  bound  set  as  per  the  following  analysis. 

Pump  Rates  for  an  Adequate  Clean  Flow  Approximation 

The  first  step  in  making  the  clean  flow  approximation  useful  is  the  determination 
of  a  rate  for  which  it  becomes  adequate.  Ultimately,  this  solution  can  always  be  used  as  a 
lower  bound  to  y  for  any  non-zero  pump  rate,  but  it’s  illustrative  usefulness  for  jc  is  not 
realized  unless  the  pump  rate  is  raised  sufficiently  high.  Figure  6.2  shows  a  pump  rate 
which  makes  the  clean  flow  approximation  of  x  poor,  while  Figure  6.3  illustrates  a  pump 
rate  yielding  a  better  approximation,  deemed  adequate.  Because  this  solution  over- 
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Clean  Flow  Approximation  (Unimodal) 
Analytic  Solution  (One  Parameter) 
Numeric  Solution  (Unimodal) 


Figure  6.2  Pump  Rate  (p=l)  Yielding  a  Poor  Clean  Flow  Approximation 


Clean  Flow  Approximation  (Unimodal) 
Analytic  Solution  (One  Parameter) 
Numeric  Solution  (Unimodal) 


t 

Figure  6.3:  Pump  Rate  (p=10)  Yielding  an  Adequate  Clean  Flow  Approximation 
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estimates  the  transfer  from  y  to  x,x  will  be  high  early  on  and  low  for  the  rest  of  the  time. 
The  criterion  for  this  determination  of  adequacy  was  that  they  were  visibly  similar 
through  the  advection  domination  stage  and  until  the  unimodal  case  diverges  from  the 
one-parameter  case. 

Comparison  of  x  and  v  for  Unimodal  Case  vs  One-parameter  Case 

Mobile  Region  Concentration.  The  fact  that  the  unimodal  case  diverges  from  the 
one-parameter  case  is  a  major  observation.  In  Figure  6.3,  the  pump  is  run  long  enough  to 
clean  up  a  major  portion  of  the  fast  and  median  sites  associated  with  the  unimodal  case 
and  a  major  portion  of  the  contaminant  in  the  one-parameter  case.  Figure  6.4  illustrates 
the  same  scenario,  but  the  pumps  ran  twice  as  long.  Pump  rate  was  10  unless  otherwise 


Clean  Flow  Approximation  (Unimodal) 
Analytic  Solution  (One  Parameter) 
Numeric  Solution  (Unimodal) 


Figure  6.4:  Longer  Run  (p=10)  For  Clean  Flow  Approximation 
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specified  in  this  entire  section.  The  divergence  begins  to  occur  at  time,  t=3,  as  slow  sites 
associated  with  the  unimodal  case  continue  their  slow  cleanup.  Note  that  for  this  given 
distribution,  nearly  three  orders  of  magnitude  difference  is  seen  between  the  two  cases  by 

t=l4r. 

The  exact  point  of  this  divergence  varies  greatly,  and  the  parameters  affecting  this 
point  are  of  great  importance.  This  point  ultimately  determines  when  consideration  of  a 
range  of  rate  parameters,  if  actually  present,  is  worth  the  effort.  If  divergence  does  not 
occur  until  safe  concentrations  are  reached  with  no  possibility  of  rebound,  then  multi¬ 
parameter  rate-limited  transport  models  are  not  worth  the  added  complexity  and 
computing  time.  On  the  other  hand,  if  this  point  is  reached  not  long  after  breakthrough  in 
a  long-term  cleanup  operation,  multi-parameter  rate-limited  transport  models  of  this  sort 
would  be  needed.  Also,  even  if  safe  concentrations  were  reached  before  this  point,  the 
possibility  of  rebound  must  be  thoroughly  investigated  before  this  sort  of  rigorous 
(although  first-order)  rate-limited  transport  modeling  is  deemed  unnecessary. 

Immobile  Region  Concentration.  Even  though  x  is  typically  what  is  predicted  by 
models  and  observed  in  the  field  and  laboratories,  the  more  important  variable  is  y.  In 
this  formulation,  it  represents  contaminant  that  is  still  stored  in  the  ground  and  in  need  of 
removal.  Only  when  contaminant  redistribution  is  examined  later,  will  y(k,t)  be 
examined  directly.  The  average  immobile  region  concentration  is  defined  as: 

KO  =  Jo  f{^)yiKt)dk .  (6.2) 

It  is  this  quantity,  when  multiplied  by  total  immobile  zone  water  content,  (0,m)totai ,  that 
quantifies  total  aqueous  contaminant  still  in  immobile  regions.  Based  upon  the  known 
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standard  deviation  of  the  Gamma  distribution,  the  code  in  Appendix  G  discretizes  X 
based  upon  a  selected  number  of  bins,  with  the  upper  bound  being  three  standard 
deviations  above  the  greater  of  the  two  means  of  the  bimodal  distribution.  The  composite 
midpoint  numerical  quadrature  approximation  to  Equation  6.2  is: 

b 

(6.3) 

«=1 

where  b  is  the  number  of  bins  and  AA,,-  is  the  bin  width.  The  following  similar  technique 
was  used  to  approximate  j(A,,t): 

y(X,t)  =  y(X,0)e~^^  +f^x(s)Xe~^^^~^^ds, 

n  (5  4) 

’  fn)  ~  X^t  hx  Xi 

m=l 

where  Wm  is  the  composite  trapezoid  weights  (lA,!,.. .,1,1/2).  In  this  work,  twenty-four  bins 
were  used. 

Figure  6.5  illustrates  the  behavior  of  y{t)  for  the  same  scenario  as  Figure  6.3. 
Note  that  the  divergence  of  the  unimodal  case  from  the  one-parameter  case  occurs  earlier 
for  y{t)  than  for  x.  Once  y(t)  begins  to  diverge,  it  takes  a  while  for  x  to  react.  Note  also 
that  the  further  out  in  time  you  go,  the  better  the  limiting  pump-rate  solution  gets. 

Change  in  Average  Effective  Rate  Parameter  with  Time 

One  very  important  process  that  is  occurring  for  the  unimodal  case  that  has  not 
been  examined  yet,  is  a  change  in  the  effective  mean  rate  parameter  with  time.  The 
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Clean  Flow  Approximation  (Unimodal),  y(t) 
Analytic  Solution  (One  Parameter),  y(t) 
Numeric  Solution  (Unimodal),  y(t) 


t 


Figure  6.5.  Behavior  of  y(t)  and  y(t)  for  One-Parameter  and  Unimodal  Cases 


expression  defining  this  parameter  is; 


W) 


(6.5) 


The  integral  in  the  numerator  was  approximated  for  the  numerical  solution  found  by  the 
code  in  Appendix  G  using  the  same  basic  technique  described  by  Equation  6.4.  The 
analytical  values  of  the  integrals  for  the  clean  flow  approximation  were  found  using 
techniques  described  in  Appendix  C.  Figure  6.6  illustrates  X(t) . 

If  Equation  6.5  could  be  used  to  obtain  a  simple  expression  for  single  rate 
parameter  variation  with  time  for  certain  general  classes  of  contaminants,  this  could  be 
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Clean  Flow  Approximation  (Unimodal) 
Numeric  Solution  (Unimodal) 


t 


Figure  6.6.  X{t)  vs  Time. 


very  useful  to  present  modelers.  Their  numerical  solution  techniques  might  allow  the 
simple  insertion  of  this  function  in  place  of  their  single  rate-parameter.  I  anticipate  that  if 
the  expression  derived  for  this  function  was  inserted  into  the  one-parameter  case 
governing  equations,  a  numerical  solution  of  them  would  closely  approximate  the  derived 
limiting  rate  solution.  A  thorough  test  of  this  hypothesis  is  warranted  before  using  this 
expression  (or  a  suitably  derived  approximation  to  it)  in  more  complex  governing 
equations.  Of  course,  the  form  of  the  theoretical  distribution  (continuous  or  discrete) 
determines  the  complexity  of  the  resulting  function.  With  this  apparent  change  in 
effective  mean  rate  parameter  in  view,  we  will  now  examine  the  soak  phase  in  detail. 
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Tailing  and  Rebound  of  Unimodal  Case  vs  One-Parameter  Case:  The  Soak  Phase 

Short  Pump,  Short  Soak.  Here,  short  is  determined  in  relation  to  the  size  of  the 
inverse  mean  rate  parameter  of  the  unchanging  distribution,  1/|X  (not  the  effective  mean 
rate  parameter  previously  described).  Figure  6.7  illustrates  a  relatively  short  pump  cycle 
with  an  equal  soaking  period.  Due  to  the  presence  of  some  faster  sites  in  the  unimodal 
case  than  the  one  parameter  case’s  sites,  the  unimodal  case  cleans  up  faster  in  this  time 


Analytic  Solution  (One  Parameter) 
Numeric  Solution  (Unimodal) 


Figure  6.7.  Behavior  of  x  for  Short  Pump,  Short  Soak. 


frame  and  appears  to  rebound  quite  similarly.  But  because  the  differences  are  quite  small, 
we  will  move  on  to  a  mid-range  pumping  period. 
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Short  Pump.  Mid-Range  Soak.  Mid-range  here  is  defined  as  approximately  ten 
inverse  mean  rates,  10/|i..  Figure  6.8  illustrates  a  parallel  sequence  of  plots  of  x,  y{t),  and 

A,(t)  for  a  one  time-unit  pump  period  followed  by  a  thirteen  time  unit  soak  period. 

Remember  that  the  analytic  solution  for  only  valid  for  the  initial  pump  period 

(with  sufficiently  high  pump  rates).  Here  we  see  x  for  the  unimodal  case  starting  below  x 
for  the  one-parameter  case,  yet  rebounding  above  it  and  continuing  to  rebound  beyond  the 
time  frame  shown.  From  Figures  6.8b  and  6.8c,  we  can  see  that  y(t)  will  clearly  drive  x 

to  a  higher  final  (asymptotic)  value  for  the  unimodal  case  and  will  take  much  longer  to 
approximately  equilibrate.  The  one-parameter  case  has  visibly  equilibrated  by  six  time 
units,  and  the  unimodal  case  is  visibly  very  close  to  it  by  fourteen  time  units. 

Even  though  x  for  these  two  cases  is  still  visibly  quite  similar.  Figure  6.8c  gives  a 
very  important  insight.  Note  that  for  over  three  time  units  after  pump  shutoff,  the 
effective  mean  rate  parameter  continues  to  decrease.  This  is  because  all  of  the  y’s  above 
X  continue  to  decrease,  and  for  this  short  time,  dominate  contaminant  redistribution  of 
contaminant  back  into  faster  sites  that  begins  at  pump  shutoff.  Even  after  fourteen  time 
units,  the  effective  rate  has  not  reattained  unity  (all  y’s  have  not  sufficiently  approached  x 
yet).  This  is  a  simple,  but  important  point  that  x  and  y  asymptotically  approach  each  other 
during  the  soak  phase.  Finally,  note  that  even  though  the  fastest  pump  rate,  if  used 
continuously,  provides  the  fastest  cleanup,  the  necessity  of  treating  renders  soaking 
beneficial  (if  the  time  can  be  afforded).  Slow  sites  continue  to  be  reduced  in 
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Clean  Flow  Approximation  (Unimodal) 
Numeric  Solution  (Unimodal) 


Figure  6.8c.  Behavior  of  A,(t)  for  Short  Pump,  Mid-Range  Soak. 


concentration  by  contaminant  movement  to  both  the  mobile  region  and  faster  sites.  We 
will  see  more  of  this  later. 

Mid-Range  Pump  and  Soak.  Figure  6.9  illustrates  this  sequence  of  a  fourteen 
time  unit  pump  and  an  equal  soak  time.  The  longer  and  harder  one  pumps,  the  larger 
these  two  cases  diverge:  over  two  orders  of  magnitude  in  this  sequence.  Note  that  the 
effective  mean  rate  gets  quite  small.  In  fact,  this  mean  rate  is  so  small  that  a  finer 
discretization  of  A,  in  the  code  would  probably  be  necessary  if  greater  accuracy  was 
desired.  Nevertheless,  a  continued  decrease  in  this  mean  rate  is  observed,  and  its  log- 
curvature  appears  to  be  decreasing.  Theory  demands  that  as  time  approaches  infinity,  this 
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Analytic  Solution  (One  Parameter) 
Numeric  Solution  (Unimodal) 


-  Analytic  Solution  (One  Parameter),  y(t) 

_ Numeric  Solution  (Unimodal),  y(0 
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Figure  6.9b.  Behavior  of  y{f)  and  y{t)  for  Mid-Range  Pump  and  Soak. 
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Figure  6.9c.  Behavior  of  A,(?)  for  Mid-Range  Pump  and  Soak. 

mean  rate  becomes  infinitely  small,  because  contaminant  would  then  only  be  present  at 
these  sites  with  an  infinitely  small  rate  parameter. 

Note  also  that  equilibration  time  remains  about  three  time  units  for  the  one-rate 
case  while  equilibration  continues  well  beyond  fourteen  soaking  time  units  for  the 
unimodal  case.  The  rebound  of  the  effective  mean  rate  is  slower  in  this  sequence  because 
slower  sites  have  now  been  somewhat  emptied,  requiring  longer  time  periods  for 
equilibration,  and  hence  longer  periods  for  rebound  of  this  rate. 

Mid-Range  Pump  and  Long  Soak.  Figure  6.10  illustrates  this  sequence  in  which 
the  soak  period  clearly  allows  x  and  y{t)  for  the  unimodal  case  to  nearly  equilibrate.  The 
apparent  slope  discontinuities  are  merely  due  to  a  numerical  output  that  was  too  coarse 
(although  step  size  remained  at  0.01).  Figure  6.10c  indeed  shows  that  our  A,(t)  has  nearly 


Clean  Flow  Approximation  (Unimodal) 
Numeric  Solution  (Unimodal) 


6.17 


Analytic  Solution  (One  Parameter) 
Numeric  Solution  (Unimodal) 


Figure  6.10a.  Behavior  of  x  for  Mid-Range  Pump,  Long  Soak. 


Analytic  Solution  (One  Parameter),  y(t) 
Numeric  Solution  (Unimodal),  y(t) 


Figure  6.10b.  Behavior  of  y(t)  and  y(t)  for  Mid-Range  Pump,  Long  Soak. 


Clean  Flow  Approximation  (Unimodal) 
Numeric  Solution  (Unimodal) 
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Figure  6. 10c.  Behavior  of  X(t)  for  Mid-range  Pump,  Long  Soak. 

reattained  unity  (the  clearest  sign  of  equilibration).  Although  longer-term  sequences  were 
investigated,  no  further  insights  were  gained. 

A  Lower  Pump  Rate  (p=l)  for  Mid-range  Pump  and  Soak.  The  sequence  shown 
in  Figure  6.11  was  produced  in  a  fashion  similar  to  Figure  6.9  to  illustrate  the  effects  of  a 
lower  pump  rate.  Note  that  the  divergence  between  the  two  cases  is  slightly  reduced 
(from  over  two  orders  of  magnitude  in  x  to  less  than  two).  Rebound  and  equilibration 
time  both  are  greatly  reduced  for  both  cases  because  a  slower  pump  rate  allows  slower 
sites  to  equilibrate  with  pumps  running,  such  that  there  is  less  equilibration  necessary 
upon  turning  the  pumps  off.  Note  that  the  limiting  rate  solution’s  effective  mean  rate  is 
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Figure  6. 1  la.  Behavior  of  x  for  p=\.  Mid-range  Pump  and  Soak. 
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Figure  6.1  lb.  Behavior  of  y{t)  and  y(t)  for  p=l,  Mid-range  Pump  and  Soak. 
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Figure  6.1  Ic.  Behavior  of  X{t)  for p=\,  Mid-range  Pump  and  Soak. 

much  lower  than  the  numerical,  further  demonstrating  that  this  pump  rate  makes  the  clean 
flow  approximation  inadequate. 

An  Approximation  to  the  Effective  Mean  Rate  Parameter.  As  a  final  note  in  this 
soak  phase  investigation,  Figure  6. 1  Ic  has  the  following  decreasing  exponential  plotted  in 
the  pumping  region: 

X{t)=ae~^\  (6.6) 

where  a  and  b  are  constants  to  be  found  based  upon  a  hypothesized  distribution  present  in 
any  given  situation.  Approximating  the  rising  rate  parameter  during  soaking  would 
require  more  investigation.  This  is  suggested  as  a  possible  rough  approximation  to  the 
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shape  of  this  effective  mean  rate’s  change  with  time.  Other  distributions  would  probably 
require  different  forms. 

Invalid  Predictions  Using  Short-Term  Data:  the  Bimodal  vs  Unimodal  Case 

In  this  section  we  will  investigate  the  effects  of  errors  in  assumed  distributions 
present  in  the  soil.  We  continue  to  include  the  one-parameter  case  also.  If,  for  example, 
the  distribution  shown  as  a  solid  line  in  Figure  6. 12  is  the  actual  distribution  present,  what 
are  the  effects  of  presuming  the  presence  of  the  other  three  distributions  shown  in  the 
figure  (each  has  a  mean  of  unity)?  Figure  6.13  illustrates  that  each  distribution  yields  a 
very  similar  solution  up  to  four  time  units.  At  this  point  the  unimodal  and  bimodal  cases 
diverge  from  the  one  rate  case.  Note  that  the  sharper  peaked  distribution  (the  lesser 
amount  of  slow  sites)  diverges  the  least  from  the  one-parameter  case.  The  unimodal 
standard  and  our  actual  distribution  continue  to  be  similar  all  the  way  until  eight  time 
units  (twice  as  long  as  it  took  to  diverge  from  the  one  rate  case). 

Let  us  assume  that  the  soil’s  rate-limited  transport  characteristics  are  to  be 
determined  in  a  laboratory  experiment  that  is  to  last  only  seven  time  units  (a  decrease  in 
concentration  at  the  well  head  by  a  factor  of  lO'"^).  Clearly,  the  presence  of  the  sharp  peak 
of  slow  sites  would  go  unnoticed  in  the  lab  and  cause  unexpected  tailing  and  rebound  in 
the  field.  The  lab  work  might  still  reveal  the  presence  of  a  distribution,  but  unless  the 
work  was  done  over  a  long  enough  period  of  time,  long  term  rate-limited  transport 
behavior  would  still  be  unpredictable  without  knowledge  of  these  slow  sites.  I  will 
suggest  an  approach  that  could  reduce  the  necessary  lab  time  later. 
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Remember  that  the  Gamma  distribution  was  chosen  for  its  flexibility  as  well  as  its 
integrability.  It  allows  for  an  assumed  non-zero,  minimum  value  for  the  rate.  Figure  6.14 
illustrates  minimum  values  of  0.25, 0.5,  and  0.75  while  keeping  a  unity  mean.  Figure 
6.15  shows  solutions  for  each  of  these  distributions  with  some  interesting  insights. 
Previously  we  noted  that  the  one  parameter  case  yielded  a  log-linear  slope  in  x  after  rate- 
limited  transport  domination.  We  also  noted  that  certain  non-zero  minimum  rate 
parameter  values  (non-zero  y‘s)  could  allow  x  to  reach  a  log-linear  slope  in  a  given  time 
period.  Figure  6.15  demonstrates  this.  Nevertheless,  note  that  the  logarithmic  slope  is 
still  less  and  concentrations  are  still  higher  when  this  log-linearity  begins.  If  there  is  a 
maximum  aggregate  size  or  a  maximum  tortuous  pore  length  into  which  contaminant  can 
diffuse,  this  would  set  a  theoretical  minimum  mass  transfer  rate.  I  therefore  suggest  using 
an  estimated  or  experimentally  derived  minimum  rate  in  this  case.  I  do  not  understand 
the  chemical  conceptualization  well  enough  to  discuss  the  merit  of  a  minimum  rate  for 
that  model. 

Redistribution  of  Contaminant  During  Soak  Phase 

In  the  soak  phase  section,  redistribution  was  mentioned  often  but  never  really 
illustrated.  Figure  6.16  actually  shows  this  redistribution.  The  horizontal  lines  are  the 
presumed  single  immobile  concentration  driven  by  a  single  rate  (a  line  instead  of  a  point 
at  X,=l  for  illustrative  purposes),  while  the  curves  represent  y  for  each  discretized  rate  and 
for  a  few  selected  times.  A  discretized  rate  here  is  referring  to  the  number  of  discrete 
values  of  X  used  to  produce  these  plots  (24).  A  fourteen  time  unit  pump  and  a  fourteen 
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Figure  6.16.  Contaminant  Redistribution. 


time  unit  soak  was  used  for  this  figure,  and  snapshots  of  y  are  shown  at  t=7, 14,21  and  28. 
Remember  that  y  is  unity  initially  for  all  rates. 

The  sharp  elbow  in  each  curve  represents  the  break  point  between  rates  for  which 
the  LEA  would  be  valid  and  those  for  which  it  would  not  be  valid.  The  sites  for  which 
the  LEA  would  not  be  valid  are  the  ones  that  cause  the  non-loglinear  (decreasing 
logarithmic)  slope  in  x.  Although  this  curvature  will  be  present  out  to  infinity  (for 
distributions  with  y=0),  it  will  become  negligible  eventually  due  to  the  small  water 
content  fractions  of  such  small  rates  (according  to  the  distribution).  The  curvature  will  be 
absent  eventually  for  non-zero  values  of  y,  as  was  previously  shown.  For  t=l,  rates  above 
approximately  1.7  appear  to  have  their  respective  y  values  equal  to  x,  while  rates  above 
1.0  (the  distribution  mean)  appear  to  have  this  condition  for  t=\A.  This  demonstrates  for 
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this  distribution  that  well  before  fast  sites  (above  the  mean)  reach  a  level  validating  the 
LEA,  unimodal  solutions  have  already  greatly  diverged  from  the  one  rate  case. 

Note  primarily  the  clear  movement  of  slow  site  contaminant  to  fast  sites  during 
the  soaking  phase  from  14  to  21  to  28.  Total  conteiminant  present  in  any  given  range  of 
sites  is  determined  by  Equation  6.2,  with  bounds  of  integration  appropriate  to  the  site 
range  of  interest.  The  observed  rates  of  redistribution  start  out  quite  fast,  but  quickly 
slow  as  contaminant  is  drawn  from  slower  and  slower  sites. 

Pulsed  Pumping  Concepts 

Now  we  seek  a  clearer  picture  of  soak  benefits.  For  this  example,  we  will 
consider  only  our  standard  unimodal  case  and  our  one  rate  case.  Figure  6.17  illustrates 
two  sequences:  a  seven  time  unit  pump  and  soak  double  cycle  versus  a  fourteen  time  unit 
pump  and  soak  single  cycle.  Note  that  total  pumping  and  soaking  time  is  identical. 

First  of  all,  a  reduced  final  concentration  is  attained  by  two  cycles  instead  of  one 
for  either  contaminant  distribution  case.  The  concentration  difference  is  most 
pronounced  for  the  one-parameter  case  because  equilibration  time  is  so  much  shorter.  In 
order  for  similar  differences  to  appear  for  the  unimodal  case,  soaking  times  must  be  long 
enough  to  allow  more  total  equilibration.  In  any  event,  a  benefit  is  still  illustrated  for  the 
unimodal  case.  Although  y  shows  very  little  difference,  jc  shows  a  significant  difference 
and  the  effective  mean  rate  plot  shows  a  marked  difference.  The  higher  effective  ending 
rate  for  the  one  cycle  case  indicates  that  less  slow  contaminant  was  cleaned  up  than  for 
the  two  cycle  case. 
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Figure  6.17c.  Behavior  of  X.(?)  for  Pulsed  Pumping. 


Finally,  although  it  is  not  easily  seen  in  Figure  6.17b,  the  two  cycle  case  has 
slightly  higher  average  concentrations  over  the  entire  pumping  regions  than  the  one  cycle 
case.  This  would  benefit  any  necessary  treating  operations. 
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VII.  Conclusions  and  Recommendations 


A  Continuous  Distribution  of  Rate  Parameters  Models  Tailing.  Rebound  and  the  Time 
and  Pump  Rate  Dependence  of  a  Single  Effective  Rate  Parameter  Well 

Rate  limited  transport  due  to  a  range  of  mass  transfer  rates,  if  present  instead  of  a 
single  rate,  is  a  very  plausible  explanation  for  many  unpredictable  behaviors.  The 
existence  of  a  distribution  of  mass  transfer  rates  would  explain  tailing  and  rebound  above 
most  rate-limited  transport  model  predictions.  More  importantly,  the  time  and  pump-rate 
dependence  of  the  single  “representative”  rate  used  by  most  present  modelers  is  implicit 
in  this  model. 

The  solutions  to  the  developed  model  clearly  demonstrate  tailing  and  rebound 
above  the  one-parameter  case,  and  a  time  and  pump  rate  dependent  effective  single  rate 
parameter.  If  rate-limited  transport  is  indeed  the  cause  of  many  instcinces  of 
unpredictable  behavior,  a  more  accurate  way  of  modeling  it  might  be  needed.  When 
more  than  one  rate-limiting  process  is  present,  this  work’s  method  is  offered  as  a  more 
realistic  (based  on  theory)  and  probably  more  accurate  means  of  modeling  rate-limited 
transport  than  most  common  methods. 

Explicit  Knowledge  of  Slow  Sites  Present  and  Deposition  History  are  Essential  to 
Accurate.  Predictive  Modeling 

Assuming  the  coexistence  of  multiple  mass-transfer  rates,  accurate,  long-term 
field  and  laboratory  predictions  require  detailed  knowledge  of  existing  slow  mass  transfer 
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rates  and  deposition  history.  It  was  clearly  seen  that  y  coupled  with  the  distribution 
determined  and  if  the  starting  values  of  y  are  not  well  known  (from  deposition  history), 
then  the  model  is  without  initial  conditions.  Further,  (even  with  accurate  initial 
conditions),  accurate  long-term  predictions  are  impossible  without  knowledge  of  the  mass 
transfer  rates  in  the  slow  site  regime  (i.e.  with  only  short  term  data). 

It  is  also  clearly  seen  in  the  literature  that  three  dimensional  models  are  inherently 
challenging.  Before  the  challenge  can  adequately  be  met,  at  the  very  least,  reasonably 
good  initial  conditions  and  good  long-term  soil  sorption/transport  characteristics  must  be 
known.  Soil  sorption/transport  characteristics  often  realistically  involve  multiple  mass 
transfer  rates.  It  is  clear  that  if  the  range  of  rates  is  wide  enough,  they  are  best  modeled 
by  either  a  discrete  or  continuous  distribution  of  rate  parameters  instead  of  just  one 
“representative”  rate  parameter.  This  work  gave  some  insight  into  how  wide  is  wide 
enough.  It  appears  that  unless  the  minimum  mass  transfer  rate  present  is  very  near  the 
representative  rate,  significant  predictive  errors  are  likely  in  the  long-term  (see  Figures 
6.14  and  6.15). 

Gaining  This  Knowledge  Probably  Requires  Long-Term  Laboratory  Experiments 

As  we  examined  the  bimodal  and  varied  unimodal  cases  for  a  hypothetical 
scenario,  it  was  clear  that  predictions  based  upon  short  term  experiments  are  only  worth 
making  for  short  term  field  situations.  Slow  sites  that  are  yet  unobserved  in  the  lab  could 
easily  be  present  in  the  field.  Also,  because  deposition  history  in  the  field  is  at  best  only 
vaguely  known,  very  good  sampling  techniques  and  tricky  laboratory  work  indeed  are 
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necessary  to  obtain  both  soil  sorption/transport  characteristics  and  initial  conditions 
simultaneously. 

One  suggested  means  of  minimizing  lab  time  necessary  to  discover  long-term  soil 
sorption/transport  characteristics  is  the  discovery  of  a  predictable  general  shape  of  the 
distribution  for  classes  of  contaminants  and  soils.  If  the  shape  was  known,  the  faster  sites 
discovered  by  short-term  experiments  could  be  used  to  estimate  the  slower  sites  using  the 
general  shape  mentioned  above.  Also,  models  which  inherently  have  the  ability  to 
estimate  parameters  from  first  principles  (such  as  diffusion-based  models),  rather  than 
repeated  curvefitting,  are  probably  closest  to  the  goal  of  accurate  predictions.  These  types 
of  models  have  the  immense  benefit  of  foregoing  long-term  laboratory  work  for  long¬ 
term  cleanups  (assuming  initial  contaminant  presence  is  known). 

Pulsed  Pumping  Generally  Appears  to  Benefit 

This  work  did  minimal  investigation  of  pulsed  pumping.  It  was  generally  found 
to  benefit  when  the  treatment  process  was  a  limiting  cost  consideration.  This  conclusion 
is  only  true  from  the  perspective  of  rate-limited  transport,  since  other  transport  processes 
were  generally  ignored.  Soak  time  must  clearly  be  available  for  pulsed  pumping  to  even 
be  considered.  In  other  words,  cleanup  will  always  go  more  slowly  using  pulsed 
pumping,  and  unless  this  longer  cleanup  time  is  acceptable,  the  efficiency  benefits  are  out 
of  reach.  This  technique  generally  caused  higher  average  concentrations  while  the  pumps 
were  running,  yielding  a  more  efficient  treatment  process.  The  efficiency  of  a  given 
treatment  process  is  maximized  normally  with  maximum  contaminant  concentrations 
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being  extracted  from  the  wells.  Also,  the  pumping  efficiency  is  increased  as  well  because 
contaminant  extracted  per  pumping  unit  time  is  increased. 

Suggestions  for  Further  Research 

First  of  all,  a  wider  variation  in  input  parameters  p  and  v  is  very  much  in  order. 
Exploring  every  conceivable  possibility  in  the  field  would  likely  give  more  insight  into 
exactly  when  such  a  model  is  necessary.  Comparing  laboratory  and  field  data  to  this 
model  could  give  insights  both  into  the  accuracy  of  its  sorption/transport  modeling  ability 
as  well  as  what  process  causes  what  effect  in  BTC’s. 

Another  beneficial  addition  to  this  work  would  be  to  allow  contaminanted  water 
to  flow  into  the  region  of  interest,  the  typical  laboratory  procedure  for  generating  BTC’s. 
This  would  make  the  above  suggestion  feasible.  Deposition  modeling  is  equally 
important. 

Actually,  this  model  without  modification  could  be  used  to  analyze  stagnant 
deposition.  A  few  computer  runs  with  this  initial  condition  were  done,  but  again,  a  much 
more  thorough  investigation  is  called  for.  As  with  cleanup,  a  thorough  investigation  of  a 
wide  range  of  possible  input  parameter  values  and  initial  conditions  would  allow  the 
application  of  many  conclusions  made  in  this  work  to  any  given  field  situation.  It  would 
also  lead  to  more  insight  into  when  such  a  model  is  indeed  necessary. 
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Appendix  A:  A  Change  of  Variables  Facilitating  Mathematica  Integrations  of  Gamma 

Distribution  Function  Moments. 


Recall  that  Equation  3.7  gives  a  standard  form  of  the  three-parameter  Gamma 
Distribution: 


1 

fk-y] 

a-1 

'-{k-y)' 

|3r(a) 

[  P  J 

exp 

L  P  J 

-oo<y<oo,  X>y,  a>0,  |3>0, 


(3.7) 


For  this  work,  y  must  always  be  greater  than  zero  (rate  constants  are  always  positive), 
which  is  why  z  is  greater  than  zero  also.  The  distribution  in  this  form  could  not  be 
integrated  by  Mathematica,  but  with  the  following  coordinate  transformation,  it  became 
integrable: 


f{z)  = 


pr(a)^p 


exp 


-(z) 

P  . 


(A.1) 


0  <  z  <  oo,  a  >  0,  p  >  0,  dz  =  dk 


The  following  Mathematica  sequence  shows  how  analytic  forms  for  p,  a,  j,  and  w  were 
obtained  (note  that  cx=a,  ^=b,  and  y=g): 


fl  [a_,b_,z_]:=(l/(b*Gamma[a]))*(z/b)^(a-l)*Exp[-z/b] 

lntegrate[fl[a,b,z],{z,0,lnfinity}]  (*Normalized  Probability  Distribution  Proven*) 
1 


A.l 


Integrate[fl [a, b,z]*(z+g),{z,0, Infinity}]  (*X.=z+g,  Proven  Distribution  Mean*) 
ab  +  g 

Sqrt[Integrate[fl[a,b,z]*(z-a*b)^2,{z,0,Infinity)]]  (*Proven  Standard  Deviation*) 
(a*b^2)^(l/2) 

Integrate[f  1  [a,b,z] *(z+g)^2*Exp[-t(z+g)] ,  { z,0, Infinity } ]  (* 

((b^(-l))^a*(a*b^2  +  a''2*b^2  +  2*a*b*g  +  g^2  +  2*a*b^2*g*t  +  2*b*g^2*t  + 
bA2*gA2*tA2))/(EA(g*t)*(b^(-l)  +  t)Aa*(l  +  b*t)^2) 

Integrate  [f  1  [a,b,z]  *  (z+g)*Exp[-t(z+g)] ,  { z,0,lnfinity }  ] 

((b^(-l))^a*(a*b  +  g  +  b*g*t))/(E''(g*t)*(b^(-l)  +  t)^a*(l  +  b*t)) 
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Appendix  B:  Mathematica  Analytic  Solution  of  One-Parameter  Case 


In  this  appendix,  the  Mathematica  sequence  used  to  obtain  the  analytic  solution  of 
the  one  rate  case  is  described.  Also  included  is  the  module  built  to  allow  pump  cycling 
with  this  solution.  Although  the  output  was  a  bit  messy  and  the  simplification  process 
was  a  little  tedious,  Mathematica  definitely  made  the  overall  task  easy: 

eqns={x'[t]+p*x[t]==v*k*(y[t]-x[t]),  y'[t]==k*(x[t]-y[t]),y[0]==y0,x[0]==x0}; 

Simplify  [DSolve[eqns,  { x,y }  ,t]] 

Without  the  simplify  command,  five  pages  of  solution  was  output.  With  it,  the  solution 
was  a  page  long,  with  the  rest  of  the  simplifying  done  by  hand.  All  hand  simplifying  was 
tested  against  plots  of  a  function  built  from  the  actual  output  of  this  sequence.  Also,  the 
solutions  calculated  by  the  code  in  Appendix  G  were  tested  against  these  forms  also.  The 
following  is  the  form  of  the  functions  used  for  this  solution: 
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xa[u_,p_,v_,xO_,yO_,t_]  := 

(l/(2*(Sqrt[(u-p)''2+u*v*(2*u+2*p+u*v)])))* 

(xO*(-u+p+u*v+(Sqrt[(u-p)^2+u*v*(2*u+2*p+u*v)]))  - 

2*yO*u*v)*Exp[-(u+p+u*v+(Sqrt[(u-p)^2+u*v*(2*u+2*p+u*v)]))*t/2]  + 

(l/(2*(Sqrt[(u-p)^2+u*v*(2*u+2*p+u*v)])))* 

(xO*(u-p-u*v+(Sqrt[(u-p)''2+u*v*(2*u+2*p+u*v)]))  + 

2*yO*u*v)*Exp[-(u+p+u*v-(Sqrt[(u-p)^2+u*v*(2*u+2*p+u*v)]))*t/2]; 

ya[u_,p_,v_,xO_,yO_,t_]  := 

(l/(2*(Sqrt[(u-p)^2+u*v*(2*u+2*p+u*v)])))* 

(-2*x0*u  +  yO*(u-p-u*v+(Sqrt[(u-p)^2+u*v*(2*u+2*p+u*v)])))* 

Exp[-(u+p+u*v+(Sqrt[(u-p)^2+u*v*(2*u+2*p+u*v)]))*t/2]  + 

(l/(2*(Sqrt[(u-p)^2+u*v*(2*u+2*p+u*v)])))* 

(2*x0*u  +  yO*(-u+p+u*v+(Sqrt[(u-p)^2+u*v*(2*u+2*p+u*v)])))* 

Exp[-(u+p+u*v-(Sqrt[(u-p)^2+u*v*(2*u+2*p+u*v)]))*t/2]; 

% 

It  was  also  necessary,  at  times,  to  use  this  function  for  a  pump  and  soak  sequence. 
This  was  made  possible  by  the  following  function,  which  essentially  reinitializes  the 
solution  at  soaking  time: 

xcyc[u_,p_,tp_,v_,xO_,yO_,tJ  :=Which[t<=tp,xa[u,p,v,xO,yO,t] , 
t>tp,xa[u,0,v,xa[u,p,v,x0,y0,tp],ya[u,p,v,x0,y0,tp],t-tp]] 
ycyc[u_,p_,tp_,v_,xO_,yO_,t_]  :=Which[t<=tp,ya[u,p,v,xO,yO,t] , 
t>tp,ya[u,0,v,xa[u,p,v,x0,y0,tp],ya[u,p,v,x0,y0,tp],t-tp]] 
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Appendix  C:  Module  for  Exponential  Integral  Form  and  Clean  Flow  Approximation 


In  the  Solution  section,  Equations  4.15  and  4.16  defined  the  task  at  hand: 


ax  1  jix  ^  ax 

^  x’"  m-l  x’^  rn-l  x^  ^ 


(4.15) 


(4.16) 


The  following  module,  ixm  was  built  to  do  Equation  4.15’s  integration  by  parts,  utilizing 
Mathematica’s  built-in  function  to  evaluate  Equation  4.16  for  the  last  integration  by  parts: 


ixm[g_,h_,av_,mv_]  :=Module[  { cum,i } , 

(*this  finds:  Integrate[Exp[a*x]/x^m,{x,g,h}]*) 

(*g  &  h  must  be  positive,  and  mv  must  be  an  integer>0*) 
cum=ExpIntegralEi  [av*h]  -ExpIntegralEi  [av*g] ; 
Which[mv>2,Do[cum=(-l/(i-l))*(Exp[av*h]/h^(i-l)- 

Exp[av*g]/g^(i- 1  ))+av*cum/(i- 1 ),  { i,2,mv }  ] , 
mv==2,cum=-Exp[av*h]/h  +  Exp[av*g]/g  +  av*cum]; 
N[cum]] 


The  following  modules  were  used  to  actually  plot  x(t),  y(t),  and  X-(t) : 
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xlim[t_,p_,v_,d  l_,a  l_,b  l_,g  I_,a2_,b2_,g2_,x0_,y0_]  :=Module[  { w  1  i,w2i,d2=  1 -d  1 } , 
w  1  i=yO*Exp[(g  1  -p)/b  1  ]  *  (a  1  *ixm[  1 , 1  +b  1  *  t,(p-g  1  )/b  1  ,a  1 + 1  ]  + 
(gl/bl)*ixm[l,l+bl*t,(p-gl)/bl,al]); 
w2i=y0*Exp[(g2-p)/b2]*(a2*ixm[l,l+b2*t,(p-g2)/b2,a2+l]  + 

(g2/b2)*ixm[  1 ,  l+b2*t,(p-g2)/b2,a2]); 

N[Exp[-p*t]*(xO  +  v*dl*wli  +  v*d2*w2i)]]; 
yave[t_,dl_,al_,bl_,gl_,a2_,b2_,g2_,y0_]:= 

dl*((bl^(-l)ral*yO)/(E''(gl*t)*(bl^(-l)  +  t)^al)+ 
(l-dl)*((b2X-l))^a2*yO)/(E^g2*t)*(b2''(-l)  +  t)^a2); 
muave[t_,a_,b_,g_,yO_]:=((b'^(-l))^a*(a*b  +  g  +  b*g*t)*yO)/ 

(EA(g=t=t)*(bA(-i)  +  t)Aa*(i  +  b*t)); 

lambar[t_,d  l_,al_,b  l_,g  I_,a2_,b2_,g2_,y0_]  := 

(dl*muave[t,al,bl,gl,yO]  +  (l-dl)*muave[t,a2,b2,g2,y0])/ 
yave[t,d  1  ,al  ,b  1  ,g  1  ,a2,b2,g2,y0] 

These  proved  quite  useful  in  the  quick  examination  of  a  multitude  of  solution  behaviors. 
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Appendix  D:  Second  Order  Convergence  and  Uncertainty  Monitoring  for  Composite 

Trapezoidal  Method 


This  appendix  will  demonstrate  the  second  order  convergence  of  this 
implementation  of  the  composite  trapezoidal  method  and  summarize  the  uncertainty 
monitoring  principles  for  this  work.  Figure  D.l  is  a  pump  and  soak  solution  (short  term) 
for  which  uncertainty  and  convergence  is  to  be  analyzed. 


Figure  D.  1 .  Short  Term  Pump  and  Soak  Solution  (/i=0.001). 


Figure  D.2  clearly  demonstrates  the  second  order  convergence  of  the  method  since  errors 
shrunk  by  a  factor  of  100  when  h  was  shrunk  by  a  factor  of  10. 


D.l 


0.0  0.5  1.0  1.5  2.0  2.5  3.0  3.5  4.0 

t 


Figure  D.2.  Change  in  Error  for  Change  in  h  from  0.01  to  0.001. 


Note  that  uncertainty  denotes  a  simple  difference  between  the  h=0.0l  and  0.001  solutions 
(the  second  common  uncertainty  monitoring  method  mentioned  in  the  Numerical 
Methods  section).  Notice  also  from  Figure  D.2  that  error  and  estimated  uncertainty 
follow  each  other.  For  a  more  important  uncertainty  reference,  %  differences  are  plotted 
in  Figure  D.3. 


Figure  D.3.  Percent  Uncertainty  and  Error  for  /i=0.01. 


From  a  percentage  reference,  there  is  actually  divergence  observed,  but  notice  how  small 
relative  error  actually  is  to  begin  with  for  ^=0.01 :  nearly  one  thousandth  of  a  percent. 
Notice  also  that  it  does  not  ever  grow  after  the  short  initial  hump. 


Appendix  E:  Third-order  Error  of  Quadratic  Interpolation 


This  appendix  will  summarize  the  proof  that  quadratic  interpolation  involves  third 
order  error  before  it  is  employed  in  any  composite  fashion.  Quadratic  interpolation  is 
accomplished  by  the  following: 

/l/2  =-^/o+|/i (E.l) 


To  start,  the  function  to  be  interpolated  must  be  expanded  in  a  Taylor’s  series  about  the 
halfway  point  for  each  point  to  be  used  in  the  approximation: 


-h 


-h 


f  ^  f  rf/  ^2/ 

Jo  -  Jl/2  + J\/2  +  ~ — Jl/2  + - — — /1/2+- 


f\  =fll2  +‘^^/l/2  +  21 


/i72  + 


-Y 

3! 


/l/2+-  , 


fl  =  /l/2  +  /i72  +  —  f{/2  +  2 ! 


v2  J 


/i72+- 


(E.2) 


Now,  if  these  are  substituted  into  our  quadratic  interpolation  formula,  E.  1  and  coeficients 
of  the  various  orders  of  derivatives  are  found,  this  yields: 


/l/2  -  /l/2  +  Y^/i72+-" 


(E.3) 


Notice  that  the  third  order  term  is  the  first  non-zero  contribution  to  errors  in  Taylor  series 
truncation,  thus  making  the  technique  third  order  in  a  single  application. 
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Appendix  F:  Second-order  Convergence  of  Linz’s  Stated  Fourth-order  Block-By-Block 

Method 

In  the  previous  appendix,  we  saw  that  Linz’s  example  block-by-block  technique 
was  clearly  only  at  most  third  order  when  applied  to  an  example  function  that  is  not 
constant.  This  author  believes  based  upon  theory  and  the  following  example  program, 
that  because  quadratic  interpolation  is  applied  to  two  integrals  per  step,  and  for  this  work 
it  would  be  applied  in  a  composite  format  (for  Z),  the  method  reduces  to  second  order  (as 
long  as  composite  integration  errors  do  not  overwhelm  interpolation  errors). 

This  concept  is  easily  proven  for  composite  integration,  as  demonstrated  by 
several  authors.  (Young  and  Gregory,  1988:373,  Hildebrand,  1987:95)  Its  basis  there 
hints  strongly  to  its  reality  for  quadratic  interpolation,  since  Newton-Cotes  type 
integration  approximation  schemes  are  interpolatory  by  nature.  But  to  further  substantiate 
the  hypothesis,  without  delving  into  a  messy  proof,  we  will  simply  code  this  method. 

The  following  code  was  used  to  generate  the  following  tables’  values. 


PROGRAM  BTest 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

DOUBLE  PRECISION  f [ALLOCATABLE] ( : ) 

INTEGER*4  m, kp , j . n, nl , n2 

CHARACTER*?  fname 

10  FORMAT  (E11.4E2, ' , ' ,E22.16E2, • , ' ,E22.16E2, ' , ' ,E22.16E2) 

20  t=0D0 

f0=lD0 

WRITE(*,*) 'Input  h,m,kp' 

READ(*,*)  h,m,kp 
c 

c  Allocate  necessary  array  sizes 

c 

ALLOCATE  (f(0:m)) 


F.l 


c  Reads  in  output  filename 
c 

WRITE(*,*)  'Output  to  file?  ( l=yes , anything  else=no) ' 

READ ( * , * )  j  f 
IF  (jf.EQ.l)  THEN 

WRITE {*,*)' Enter  7-char  output  f ile ( [appos] f ilenam[appos] ) 
READ ( * , * )  f name 
c 

c  File  for  data  to  be  written  to 

c 

OPEN  ( 1 , FILE=  f name / / ' . out ' , STATUS= ' UNKNOWN ’ , 

+  ACCESS= ’ SEQUENTIAL ' , FORM= ' FORMATTED ' ) 

WRITE  (1,’^)  'fO,h=' 

WRITE  (1,*)  fO,h 
WRITE  (1,*) 

WRITE ( 1 , * ) ' t  ' , f name , ' . num  ' , f name , ' . ana  ' 

WRITE  (1,10)  t,f0,f0,t 
END  IF 
WRITE  ( * , * ) 

WRITE  (*,10)  t,f0,f0,t 
c 

c  Initializing  variables 

c 

f (0)=f0 
c 

c  Start  loop 

c 

100  DO  1000  n=0,m-2,2 
nl=n+l 
n2=n+2 
c 

c  Step  times 

c 

t=n*h 

th=:{n  +  0.5D0)*h 
tl=nl*h 
t2=n2*h 
c 

c  Summations 

c 

sl=0D0 

IF  (n.NE.O)  THEN 
DO  200  j=0,n 

IF  { (j .EQ.O) .OR. (j .EQ.n) )  THEN 
wt=lD0 

ELSE  IF  {MOD(j ,2) .EQ.O)  THEN 
wt=2D0 

ELSE 

Wt=4D0 
END  IF 

sl=sl  +  Wt*f ( j ) 

200  CONTINUE 

END  IF 

f(nl)=((f(n)  +  (h/6D0)*(-15D0*f (n)-20D0*sl*{h/3D0)-32D0* 

+  ((h/3D0)*sl  +  h*f (n) *5D0/12D0)+4D0* {h/3D0) * (si  +  f{n))))  + 

+  ( (h/6D0) * (3D0  -  32D0* (-h/12D0)  +  4D0*h/3D0))* 

+  (f (n) +  (h/3D0) * {-6D0*f (n) -8D0*sl* {h/3D0 ) -32D0 * 

+  ({h/3D0)*sl  +  h*f (n)*5D0/12D0)-8D0*(h/3D0)*(sl  +  f(n))))/ 

+  (lD0-(h/3D0)*{-32D0*{-h/12D0)  -  6D0  -  8D0*h/3D0 ) ) ) / 

+  (IDO  -  (h/6D0) * (-24D0  -  24D0*2D0*h/3D0  +  4D0*4D0*h/3D0 )  - 
+  ( (h/6D0) * (3D0  -  32D0* (-h/12D0)  +  4D0*h/3D0))* 

+  ((h/3D0)*(  -24D0  -  32D0*2D0*h/3D0  -  8D0*4D0*h/3D0  ))/ 

+  (IDO- (h/3D0) * {-32D0* (-h/12D0)  -  6D0  -  8D0*h/3D0) ) ) 

f (n2)=( (f (n)+  (h/3D0) * ( -6D0*f (n) -8D0*sl* (h/3D0 ) -32D0* 
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+  {(h/3D0)*sl  +  h*f {n)*5D0/12D0)-8D0*(h/3D0)*(sl  +  f{n))))  + 

+  {(h/3D0)*(  -24D0  -  32D0*2D0*h/3D0  -  8D0 *4D0*h/3D0  ))*f(nl))/ 

+  {lD0-(h/3D0)*(-32D0*{-h/12D0)  -  6D0  -  8D0*h/3D0) ) 

c 

c  Calculate  fan 
c 

fan2=2D0*DEXP(-4D0*t2)  -  DEXP ( -2D0*t2 ) 
diff2=fan2-f (n2) 
c 

c  Output  to  screen  and  file  values  for  f  each  kp_th  step 
c 

IF  (kp.EQ.l)  THEN 

fanl=2D0*DEXP(-4D0*tl)  -  DEXP ( -2D0*tl ) 
dif fl=fanl-f (nl) 

WRITE  (*,10)  tl,f (nl) ,fanl,diffl 
WRITE  (*,10)  t2,f (n2) ,fan2,diff2 
IF  (jf.EQ.l)  THEN 

WRITE  (1,10)  tl, f (nl) , fanl,diffl 
WRITE  (1,10)  t2,f (n2) ,fan2,diff2 
END  IF 

ELSE  IF  (MOD(n2,kp) .EQ.O)  THEN 

WRITE  (*,10)  t2,f (n2) ,fan2,diff2 
IF  (jf.EQ.l)  THEN 

WRITE  (1,10)  t2, f (n2) ,fan2,diff2 
END  IF 
END  IF 

1000  CONTINUE 

IF  (MOD(n2,kp) .NE.O)  THEN 

WRITE  (*,10)  t2,f (n2) ,fan2,diff2 
IF  (jf.EQ.l)  THEN 

WRITE  (1,10)  t2,f (n2) ,fan2,diff2 
END  IF 
END  IF 

DEALLOCATE  (f) 

C 

c  Allow  another  run 
c 

WRITE  (*,*)  'Another  run?  (l=yes , anything  else=no) ' 

READ  ( * , * )  mo 

IF  (mo.EQ.l)  GOTO  20 

END 


Table  F.  1.  Error  for  Block-by-block  Technique,  Non-constant  Example  (function 
_ _ location:  Linz,  1985:187) _ 


h= 

0.1 

0.0  i 

0.001 

t=0.1 

2.03E-3 

2.05E-5 

2.05E-7 

t=0.2 

1.98E-3 

2.25E-5 

2.25E-7 

This  technique  is  clearly  proven  to  be  only  second  order  for  this  example. 
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Appendix  G:  Composite  Trapezoidal  Code  for  Integro-differential  Equation,  With 
Bimodal  Case  Allowed,  Three  Different  Pump  Changing  Schemes,  Flexible  Output  ofX 
and  Y,  and  Calculation  of  Mean  Y  and  Effective  Mean  X 


c  This  program  solves  for  the  mobile  water  concentration,  X,  in  a 

c  fixed  volume,  V,  of  constant  saturation  soil.  A  single  volumetric  flow 

c  parameter,  Q,  is  used  to  extract  solvent  that  has  solute (s) 

c  dissolved  in  it.  Clean  solvent  is  assumed  to  be  flowing  through  the 

c  outer  boundaries  of  V  as  the  pump  runs.  The  soil  is  modeled 

c  as  a  continuum  of  immobile  regions  through  which  no  solvent  is  flowing, 

c  No  solute (s)  are  assumed  sorbed  within  these  regions,  and  solutes  have 

c  only  a  diffusive  path  to  follow  in  order  to  reach  the  mobile  zone, 

c  Whether  due  to  different  sorption  site  characteristics  or  different 

c  geometries  of  immobile  regions,  solute  moves  from  immobile  to  mobile 

c  regions  based  upon  a  distribution  of  rate  parameters,  lam  (lambda), 

c  having  a  probability  density  function,  f{lam).  Two  cases  are  allowed 

c  for  (initially) :  2  average  lam's,  ul  and  u2  (with  dl%  in  ul  and 

c  d2%  in  u2,  and  a  corresponding  average  immobile  concentration,  Y) , 

c  and  a  bimodal  Gamma  distribution  of  lam's  (with  dl%  in  fl(lam)  and  d2%  in 

c  f2(lam),  and  a  corresponding  continuum  of  immobile  regions  with 

c  concentrations,  Y(lam,t).  A  first  order  rate  model  is  used  for  the 

c  non- equilibrium  transfer  of  solute  to  the  mobile  region,  and  X  and  Y 

c  both  are  modeled  as  volume-averaged  concentrations.  A  ratio  of  immobile 
c  to  mobile  solvent  content  in  the  soil,  v,  is  also  used, 

c 

c  The  governing  two  differential  equations  (with  these  assumptions)  are: 
c 

c  1.  Ydot=lam* [X ( t) -Y (lam, t ) ] 
c 

c  2.  Xdot+p*X=v*Integral  from  0->infinity [f (lam) *lam* [Y(lam, t) -X(t) ] ] dlam 

c 

c  P=Q/ (V*kr) ,  where  kr  is  a  reference  rate  constant  used  to 

c  non-dimens ionalize  everything 

c  Dot  implies  the  partial  derivative  with  respect  to  time 

c 

c  The  first  equation  is  solved  by  an  integrating  factor  to  obtain: 

c 

c  Y ( lam, t) =Y ( lam, 0 ) *Exp [-k*t] + 

c  Integral  from  0->t [lam*X(t " ) *Exp [ -lam* (t-t" ) ] dt " 

c 

c  This  is  plugged  into  equation  #2  and  manipulated  algebraically  to 

c  obtain  the  form  which  is  solved  numerically, 

c 

c  Xdot  is  dealt  with  by  integration  then  trapezoid  rule,  and 
c  the  integral  from  0->t  is  dealt  with  using  composite  trapezoid  rule, 

c 

c  The  Gamma  distributions  used  have  parameters  alpha,  beta,  and  gamma 

c  (a,b,g).  Their  form  is: 

c 

c  f  (lam)  =  (l/(b*GAMMA[a]  ))*(  (  ( lam-g) /b) (a-1 )  )  *Exp  [  -  ( lam-g) /b] 

c  Integrals  involving  this  were  analytically  found  and  entered  into  the 

c  code  in  appropriate  places.  u=mean  of  each  distribution=a*b  +  g,  and 

c  the  standard  deviation  of  each  distribution  is  b*Sqrt(a) 

c 

PROGRAM  TRAPS 
C 

c  Variable  Definitions 

c 

c  General : 

c  to=previous  time 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


t=present  time 

m=total  #  of  time  steps  to  be  taken 
h=time  step  size 
no=previous  time  step 
n=present  time  step 

q=p+v*u  (P/U  defined  above, pO“>qo  from  previous  time  step) 
lam( j ) ,y { j ) =j  bins  of  discretized  lambda ' s/y ’ s , j  =  1 ... #  of  bins 
yO=starting  concentration  in  each  immobile  region 
yOa=present  yO  for  analytic  solution 
x=concentration  in  mobile  region 
x{l, i) =trapezoid  @  ith  step 

x{2 , i) =trapezoid  @  ith  step  using  average  model 

w=function  of  t  (different  when  average  lam  is  used,  called  wa) 

wf (t) =v* Integral  from  0->Infinity 

[f (lam) *lam*yO (lam) *Exp[-lam*t] ]dlam 
wfa ( t) =v*u*yO*Exp [-u*t] 

wo=wf (previous  time) .. .woa=wf a (previous  time) 
tau=t-t ' 

k(i) =kf (tau) ,  function  of  tau  (different  for  average  lam,  ka(i)) 
kf  (tau)  =v* Integral  from  0->inf  inity  [ f  (lam)  *lam'^2 *Exp  [ -lam*tau]  ]  dlam 
kfa ( tau) =v*u^2  *  Exp [-u* tau] 
xO=starting  value  of  x 

xOa=present  initial  value  of  x  for  analytic  solution 

npx=number  of  steps  between  output  of  x 

npy=number  of  steps  between  output  of  y 

nps=number  of  steps  taken  at  last  pump  rate  change 

npc=number  of  steps  to  be  taken  before  next  pump  change 

nc=number  of  pump  and  soak  cycles  to  be  taken 

jox,joy=#  of  outputs  of  x  and  y  per  pump  on  cycle 

jfx,jfy=#  of  outputs  of  X  and  y  per  pump  off  cycle 

cl-4f (p) , cmlf , cm2f=constants  for  analytic  solution  evaluation 

cl-4n, cml , cm2=present  values  for  constants  in  analytic  solution 

Minimal  variable  declarations  (remember:  allocatable  isn't  F-77) 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

x[ ALLOCATABLE] ( : , : ) 
k [ALLOCATABLE] ( : ) 
ka [ALLOCATABLE] ( : ) 
kf , kf a , lamb , lamv 
y[ ALLOCATABLE] ( : ) 
lam [ALLOCATABLE] ( : ) 
f [ALLOCATABLE] ( : ) 

m, j , ns , n , npc , npo , npf , nc , garni , gam2 , gamv 
jox, joy, j  fx, jfy,npx,npy 
nal , na2 , nav 
fname 

Functions 


DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
INTEGER *4 
INTEGER *4 
INTEGER*! 
CHARACTER* 7 


wf (tv,uv, vv,nav,bv,gv,yOv) = (uv  +  bv*gv*tv) *yOv*vv/ 

+  (DEXP(gv*tv) * ( (IDO  +  bv*tv) ** (nav  +  IDO) ) ) 

wfa ( tv, uv, w, yOv) =uv*vv*yOv*DEXP ( -uv*tv) 
kf (tau,uv, vv,nav,bv, gv) =vv* ( 

+  (uv**2D0)  -f  (nav*  (bv**2D0)  *  (IDO  +  2D0*gv*tau)  )  + 

+  ( (gv**2D0) * (tau**2D0) * (2D0*bv  +  (bv**2D0 ) ) ) ) *DEXP ( -gv* tau) / 

+  ((IDO  +  bv*tau) ** (nav  +  2D0)) 

kfa ( tau,uv, w) = (uv**2D0) *vv*DEXP ( -uv*tau) 

zf (pv,uv, w) =SQRT ( (uv-pv) **2D0  +  uv*vv* (2D0* (uv+pv)  +  uv*vv) ) 
cmlf (pv, uv, w) = (-pv  -  uv  -  zf(pv,uv,vv)  -  uv*vv) /2D0 
cm2f (pv,uv, w) = (-pv  -  uv  +  zf(pv,uv,vv)  -  uv*vv) /2D0 
clf  (pv, uv,  w, xOv, yOv)  =  (IDO/  (2D0*zf  (pv, uv,  w)  )  )  * 

+  (xOv* (pv  -  uv  +  zf(pv,uv,vv)  +  uv*vv)  -  2DO*yOv*uv*vv) 
c2f (pv,uv, vv,x0v,y0v) = (IDO/ (2D0*zf (pv,uv, w) ) ) * 

+  (x0v*(-pv  +  uv  +  zf(pv,uv,vv)  -  uv*vv)  +  2D0*y0v*uv*vv) 
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c3f  {pv,uv,  w,xOv,yOv)  =  (IDO/  {2D0*zf  (pv,uv,  w)  )  )  * 

+  (yOv*(~pv  +  uv  +  zf(pv,uv,vv)  -  uv*w)  -  2D0*x0v*uv) 
c4f  {pv,uv,  w,x0v,y0v)  =  (IDO  /  {2D0*zf  (pV;UV,  w)  )  )  * 

+  (yOv*  (pv  -  uv  +  zf(pv,uv,w)  +  uv*w)  +  2D0*x0v*uv) 
gamf  (nav,bv,  gv,  gamv,  laitiv)  =  (IDO/  (gamv*bv) )  * 

+  ( { (lamv-gv) /bv) ** (nav-lDO) ) *DEXP( (gv-lamv) /bv) 

c 

c  Formats 

c 

c  File  x's 

10  FORMAT ( ' , ' ,E11.4E2, ' , ’ ,E22.16E2, ' , ' ,E22.16E2, ' , ' ,E22.16E2) 
c  Pump  on  X ' s 

12  FORMAT  ( 'pon, ' ,E11.4E2, ' , ' ,E22.16E2, ' , ' ,E22.16E2, ’ , ' ,E22.16E2) 
c  Pump  off  x's 

14  FORMAT  ( 'poff , ' ,E11.4E2, ' , ’ ,E22.16E2, ' , ’ ,E22.16E2, ' , ' ,E22.16E2) 

15  FORMAT  (ElO . 4E2 , E13 . 7E2 , ElO . 4E2 , E13 . 7E2 , ElO . 4E2 , ElO . 4E2 , ElO . 4E2 ) 

c  Input,  header  info 

16  FORMAT  ( • ,  ' , E8 . 3E2 ,  ' ,  ' , E8 . 3E2 ,  ' ,  ' , E8 . 3E2 ,  ' ,  ' , 13 ,  ' ,  '  , 13 , 

+  ' , ' ,E8.3E2, ■ , ’ ,E8.3E2, ' , ' ,E8.3E2) 

c  File  lambda, y  format 

17  FORMAT  (E10.4E2, ' , ' ,E10.4E2) 
c 

c  Input  initial  pump  rate, pi,  distribution  parameters  for  each  of  two 

c  distributions,  dl  (%  in  1st  dist.),  nal,  na2 ,  ul,  u2 ,  (bl,  b2  calc.) 

c  gl,  g2,  initial  conditions  xO,  yO,  step  size,  h,  #  of  pump  cycles,  nc, 

c  #  of  pump  steps,  npo,  #  of  printouts  of  x  and  y  per  pump  cycle, 

c  iox,  ioy,  #  of  soak  steps,  npf,  #  of  printouts  of  x  and  y  per  soak 

c  cycle,  ifx,  ify. 

c 

20  CONTINUE 

WRITE {*,*) 'Default  distribution? { (dl , a, u, g) =(. 5 , 3 , 1 , 0) ) (1-yes) ' 

READ  (*,*)  ndef 
IF  (ndef.EQ.l)  THEN 
dl=lD0 
nal=3 
na2=nal 
ul=lD0 
u2=ul 
gi=0D0 
g2=gl 

ELSE 

WRITE  (*,*)  ’ Input  dl, nal, na2,ul,u2,gl,g2 ' 

READ  (*,*)  dl,nal,na2 ,ul, u2 , gl, g2 
END  IF 
d2=lD0-dl 
u=dl*ul  +  d2*u2 

25  WRITE  (*,*)  'Input  pi , v, xO , yO , h, j ox, joy , j fx, j fy , la ' 

READ  (*,*)  pi, v,x0,y0,h, jox, joy, jfx, jfy, la 

WRITE  (*,*)  'Set  cycle  times (=0),  thresholds ( -1 ) ,  or  %chng(=2)?' 

READ  (  * , * )  ncyc 
IF  (ncyc.EQ.O)  THEN 

WRITE  (*,*)  'Input  nc, npo, npf 
READ  (*,*)  nc, npo, npf 
m= (npo+npf ) *nc 
npc=npo 

ELSE  IF  (ncyc.EQ.l)  THEN 

WRITE  (*,*)  'Input  m,pft,pot' 

READ  (*,*)  m,pft,pot 
npc  =  0 

ELSE 

WRITE  (*,*)  'Input  m,po' 

READ  (*,*)  m,po 
npc-0 
END  IF 
P=pi 
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cln=clf (p, u, V, xO ; yO) 
c2n=c2f (p, u, V, xO , yO) 
tb= (LOG (c In) “LOG (c2n) ) /zf (p, u, v) 

WRITE  (*,*)  restart?  (l=yes) ' 

READ  (*,*)  St 
IF  (st.EQ.l)  GOTO  20 
c3n=c3f (p, u, V, xO , yO) 
c4n=c4f {p,u, v,xO,yO) 
cml=cmlf (p, u, v) 
cm2=cm2f (p, u, v) 
npx=jox 
npy=joy 
c 

c  Initial  conditions  applied  to  anaytic  solution 

c 

y0a=y0 

yd=y0 

x0a=x0 

c 

c  Initial  time  set  to  0 

c 

t=0D0 

ta=t 

c 

c  Allocate  necessary  array  sizes 

c 

ALLOCATE  (x{l:2,0:m) ) 

ALLOCATE  (k(0:m)) 

ALLOCATE  (ka(0:m)) 

ALLOCATE  (y(l:la)) 

ALLOCATE  ( lam ( 1 : la ) ) 

ALLOCATE  (f (1:1a) ) 

c 

c  Reads  in  output  filename 

c 

diff=0D0 

WRITE(*,*)  'Output  to  file?  ( l=yes , anything  else=no) ' 

READ ( * , * )  j  f 
IF  (jf.EQ.l)  THEN 

WRITE (*,*)' Enter  7-char  output  f ile ( [appos] f ilenam[appos] ) : ' 

READ ( * , * )  f name 
c 

c  File  for  x  data  to  be  written  to  with  header  info  and  initial  conditions 

c 

OPEN  ( 1 , FILE=fname/ / ' . exx ’ , STATUS= ' UNKNOWN ' , 

+  ACCESS= ' SEQUENTIAL ' , FORM= ' FORMATTED ' ) 

WRITE  (1,*)  'pi,v,dl,nal,na2,ul,u2,gl,g2,x0,y0,h=' 

WRITE  (1,*)  pi,v,dl,nal,na2,ul,u2,gl,g2;X0,y0,h 
WRITE  (1,*) 

WRITE ( 1 / * ) ' t*Kr  ' , fname ; ' . dis  ' , fname, ' . err  ' , fname , ' . ave  = ' 

WRITE  (1,10)  t,x0,diff,x0 
c 

c  File  for  y  data  to  be  written  to  with  header  info 

c 

OPEN  ( 2 , FILE=  fname / / ' . why ’ , STATUS= ' UNKNOWN ' , 

+  ACCESS= ■ SEQUENTIAL ' , FORM= ' FORMATTED ’ ) 

WRITE  (2,*)  ' pi , V, dl , nal , na2 , ul , u2 , gl , g2 , xO , yO , h= ' 

WRITE  (2,*)  pi,V;dl,nal,na2,ul,u2,gl,g2,xO,yO,h 
WRITE  {2 ,*)  ' f name= ' , fname 

WRITE  (2,*) 

END  IF 
c 

c  Header  info  to  screen 

c 

WRITE  ( * , * ) 
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WRITE  (*,*)  fnaitie,  '  ,  t,xdis,  err,xan,ydis,yan,  lanib=  ' 

WRITE  (*,15)  t,xO,diff ,xO,yd,yO,u 
c 

c  Initializing  several  variables 

c 

bl= (ul-gl ) /nal 
b2= (u2-g2 ) /na2 
q=p+v*u 

k(0) =dl*kf (t,ul,v,nal,bl, gl)  +  d2*kf (t,u2,v,na2,b2,g2) 
ka ( 0 ) =kfa ( t , u, v) 

X ( 1 , 0 ) =xO 
X  (2 , 0)  =x0 
xan=x0 

CALL  gam  (nal, garni) 

CALL  gam  (na2,gam2) 

bs=( (u2+3D0*b2*SQRT(na2) ) ~gl) /la 

DO  50  j=l,la 

lam(j)=gl  +  (2*j  -  l)*bs/2 
f (j ) =dl*gamf (nal, bl,gl, garni, lam(j) )  + 

+  d2*gamf (na2 ,b2 , g2 , gam2 , lam( j ) ) 

50  CONTINUE 
c 

c  Present  values  for  w,  wa,  and  constants  for  analytic  solution 

c 

w=dl*wf ( t , ul , V, nal , bl , gl , yO )  +  d2  *wf ( t , u2 , v, na2 , b2 , g2 , yO ) 
wa=wf a ( t , u , V , y 0 ) 
ns=0 
c 

c  This  ensures  no  integral  contributions /summations  on  the  first  step 
c 

jd=0 

c 

c  Constants  to  minimize  calculations  in  the  loop 

c 

cl=h/2D0 

c2=cl**2 

c 

c  start  loop 

c 

DO  1000  n=l,m 
c 

c  Step  time 

c 

no=n-l 

t=n*h 

c 

c  Sets  previous  "new"  values  for  w  equal  to  the  one  time  step  old  ones 

c 

wo=w 

woa=wa 

c 

c  New  values  for  w  and  wa 

c 

w=dl*wf (t,ul, v,nal,bl,gl,y0)  +  d2*wf (t,u2, V,na2,b2,g2,y0) 
wa=wf a ( t , u , V , y 0 ) 
c 

c  Present  values  for  k  and  ka 

c 

k(n) =dl*kf (t,ul, v,nal,bl, gl)  +  d2 *kf ( t , u2 , v, na2 , b2 , g2 ) 
ka (n) =kfa ( t, u, v) 
c 

c  Summations 

c 

stf=0D0 
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stb=ODO 

stfa=ODO 

stba=ODO 

IF  {n.GT.2)  THEN 
DO  200  j=l,n-2 

stf=stf  +  X ( 1 , j ) *k (no- j ) 
stb=stb  +  X (1, j ) *k (n-j ) 
stfa=stfa  +  X (2 , j ) *ka (no- j ) 
stba=stba  +  x (2 , j ) *ka (n- j ) 

200  CONTINUE 

END  IF 

stb=stb  +  jd*x (l,no) *k(l) 
stba=stba  +  jd*x(2 , no) *ka (1) 
c 

c  Trapezoid  steps 
c 

X { 1 , n) = (x ( 1, no)  +  cl*(w  +  wo  -  q*x(l,no)  + 

4-  jd*cl*  (x(l,  0)  *k  (no)  +  x(l,no)*k(0)  +  2D0*stf)  + 

+  cl*(x(l,0)*k(n)  +2D0*stb)))/ 

+  (IDO  +  cl*q  -  k(0)*c2) 

X (2 , n) = (x (2 , no)  +  cl*(wa+  woa-  q*x{2,no)  + 

+  jd*cl*  (x  (2 , 0 )  *ka  (no)  +  x  (2  ,  no)  *ka  (0)  +  2D0*stfa)-f 

+  cl* (x (2 , 0) *ka (n) +  2D0*stba) ) ) / 

+  (IDO  +  cl*q  -  ka(0)*c2) 

c 

c  First  step  complete 

c 

jd=l 

c 

c  Calculate  xan 

c 

ta= (n-ns) *h 

xan=cln*DEXP(cml*ta)  +  c2n*DEXP (cm2 *ta) 
c 

c  Pump  changing  scheme 

c 

IF  (ncyc.EQ.l)  THEN 

IF  { (p . EQ . ODO ) . AND. (x ( 1 , n) . GT .pot ) )  npc=n-ns 
IF  ( (p.EQ.pi) .AND. (x{l,n) .LT.pft) )  npc=n“ns 
ELSE  IF  (ncyc.EQ.2)  THEN 

pd=100D0*ABS (x ( 1 , no)  -  x(l/n))/xan 
IF  (pd.LT.po)  npc=n-ns 
END  IF 

IF  ( (n-ns) .EQ.npc)  THEN 
di f  f =xan-x ( 2 , n ) 

y0b=c3n*DEXP (cml*ta)  +  c4n*DEXP {cm2*ta) 
c  Calculate  y(lam,t)  for  output  time  and  specified  discrete  lam's 

yd=0 
temp=0 

DO  400  j=l,la 
sy=0 

DO  300  jl=l,n-l 

sy=sy  +  x(l, jl) *DEXP(h*lam( j ) * ( jl-n) ) 

300  CONTINUE 

y(j)={y0  +  lam(j ) *h*x0/2) *DEXP(-lam(j ) *t)  + 

+  lam(j)*h*(sy  +  x(l,n)/2D0) 

yd=yd  +  f ( j ) *y ( j ) *bs 
temp=temp  +  f ( j ) *y ( j ) *bs*lam( j ) 

400  CONTINUE 

lamb=temp/yd 

npc=0 

ns=n 

c  Reset  analytic  solution  IC's 

y0a=c3n*DEXP(cml*ta)  +  c4n*DEXP (cm2 *ta) 
x0a=xan 
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c  Reset  jp,  npc,  and  nip  is  set  to  ensure  proper  output;  also, 

c  output  has  pump  on  or  pump  off  anotation  next  to  it. 

IF  (p.EQ.ODO)  THEN 
p=pi 

IF  (ncyc.EQ.O)  npc=npo 

npx=jox 

npy=joy 

WRITE  (*,15)  t,x(l,n) , dif f ,xan,yd,yOb, lamb 
IF  (jf.EQ.l)  WRITE  (1,12)  t , x ( 1 , n) , xan, dif f 

ELSE 

p=0D0 

IF  (ncyc.EQ.O)  npc=npf 
npx= j  fx 
npy= j  fy 

WRITE  (*,15)  t , X ( 1 , n) , dif f , xan, yd, yOb, lamb 
IF  (jf.EQ.l)  WRITE  (1,14)  t , x ( 1 , n) , xan, dif f 
END  IF 

c  Reset  p,  q  and  analytic  solution  constants  for  pump  change 
cln=clf (p, u, V, x0a,y0a) 
c2n=c2f (p, u, V, x0a,y0a) 
c3n=c3f (p,u, v,x0a,y0a) 
c4n=c4f (p, u, V, x0a,y0a) 
cml=cmlf (p,u, v) 
cm2=cm2f (p,u, v) 

c  Time  shift  for  analytic  solution 
q=p+v*u 
END  IF 

IF  ( (MOD(n"ns,npx) .EQ.O) .AND. (ns.NE.n) )  THEN 
dif f=xan-x (2 , n) 

y0b=:c3n*DEXP  (cml*ta)  +  c4n*DEXP  (cm2 *ta) 
c  Calculate  y(lam,t)  for  output  time  and  specified  discrete  lam's 

yd=0 
temp=0 

DO  460  j=l,la 
sy=0 

DO  430  jl=l,n-l 

sy=sy  +  x(l, jl) *DEXP(h*lam( j ) * ( jl-n) ) 

430  CONTINUE 

y(j)=(y0  +  lam(j)*h*xO/2)*DEXP(-lam(j)*t)  + 

+  lam(j)*h*(sy  +  x(l,n)/2D0) 

yd=yd  +  f(j)*y(j)*bs 
temp=temp  +  f ( j ) *y ( j ) *bs*lam( j ) 

460  CONTINUE 

lamb=temp/yd 

WRITE  (*,15)  t,x(l,n)  ,diff  ,xan,yd,yOb,lainb 
IF  (jf.EQ.l)  WRITE  (1,10)  t , X ( 1 , n) , xan , dif f 
END  IF 

IF  (MOD (n-ns,npy) .EQ.O)  THEN 
c  Output  y 

IF  (jf.EQ.l)  THEN 
WRITE  ( 2 , * ) 

WRITE  (2,*)  'lam(j) ,y_t=' ,t 
DO  500  j=l,la 

WRITE  (2,17)  lam(j),y(j) 

500  CONTINUE 

END  IF 

IF  (jf.EQ.l)  THEN 

WRITE  (2 , *) 'y Ob (analytic) ,yd, lamb= ' ,y0b,yd, lamb 
WRITE  ( 2 , * ) 

END  IF 
END  IF 

IF  ( n . EQ . ns )  WRITE  ( * , * )  ' Pump  changed ! ' 

1000  CONTINUE 
c 

c  Deallocation 
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c 

DEALLOCATE  (x) 

DEALLOCATE  (k) 

DEALLOCATE  (ka) 

DEALLOCATE  (y) 

DEALLOCATE  (lam) 

DEALLOCATE  ( f ) 

C 

c  Allow  another  run 
c 

IF  (jf.EQ.l)  THEN 
CLOSE  (1) 

CLOSE  (2) 

END  IF 

WRITE  {*,*)  'Another  run?  { l=:yes , anything  else=no)  ' 
READ  ( * , * )  mo 
IF  (mo.EQ.l)  GOTO  20 
END 

SUBROUTINE  gam  (naif, nans) 

INTEGER* 4  nans 

INTEGER*!  nalf,m 

nans=l 

DO  2000  m=2,nalf-l 
nans=nans*m 
2000  CONTINUE 

END 


G.8 


Bibliography 


Adams,  Thomas  A.  and  Robert  C.  Viramontes.  Analytical  Modeling  of  Aquifer 

Decontamination  by  Pulsed  Pumping  When  Contaminant  Transport  is  Affected  by 
Rate-Limited  Sorption  and  Desorption.  MS  Thesis,  AFIT/GEE/ENG/93-S1. 
School  of  Engineering,  Air  Force  Institute  of  Technology  (AU),  Wright-Patterson 
AFB,  OH,  September  1993  (AD-A271  105). 

Ball,  William  P.  and  Paul  V.  Roberts.  “Long-Term  Sorption  of  Halogenated  Organic 
Chemicals  by  Aquifer  Material.  2.  Intraparticle  Diffusion,”  Enyironmental 
Science  and  Technology:  25:  1237-1249(1991). 

Brusseau,  Mark  L.  and  P.S.C.  Rao.  “Sorption  Nonideality  During  Organic  Contaminant 
Transport  in  Porous  Media,”  Critical  Reyiews  in  Enyironmental  Control.  19:  33- 
99  (1989). 

Burden,  Richard  L.  and  J.  Douglas  Faires.  Numerical  Analysis  (Third  Edition).  Boston: 
Prindle,  Weber  and  Schmidt  Publishers,  1985. 

Claborn,  Billy  J.  and  Ken  A.  Rainwater.  “Well-Field  Management  for  Energy 

Efficiency,”  Journal  of  Hydraulic  Engineering.  117:  1290-1303  (October  1991). 

Connaughton,  Doreen  F.  et  al.  “Description  of  Time- Varying  Desorption  Kinetics: 

Release  of  Naphthalene  From  Contaminated  Soils,”  Enyironmental  Science  and 
Technology.  27:  2397-2403  (1993). 

Cooney,  D.  O.  et  al.  “Effects  of  Particle  Size  Distribution  on  Adsorption  Kinetics  in 
Stirred  Batch  Systems,”  Chemical  Engineering  Science.  38:  1535  (1983). 

Cussler,  E.  L.  Diffusion:  Mass  Transfer  in  Fluid  Systems.  New  York:  Cambridge 
Uniyersity  Press,  1984. 

Domenico,  Patrick  A.  and  Franklin  W.  Schwartz.  Physical  and  Chemical  Hydrogeology. 
New  York:  John  Wiley  and  Sons,  Inc.,  1990. 


Gilbert,  Richard  O.  Statistical  Methods  for  Enyironmental  Pollution  Monitoring.  New 
York:  Van  Nostrand  Reinhold,  1987. 


Hartman,  Richard  T.  Optimal  Pulsed  Pumping  for  Aquifer  Remediation  When 

Contaminant  Transport  is  Affected  bv  Rate-Limited  Sorption:  A  Calculus  of 
Variation  Approach.  MS  Thesis,  AFrr/GEE/ENC/94S-2.  School  of  Engineering, 
Air  Force  Institute  of  Technology  (AU),  Wright-Patterson  AFB,  OH,  September 
1994  (AD-A284  700). 

Heyse,  Edward.  Mass  Transfer  Between  Organic  and  Aqueous  Phases:  Investigations 
Using  A  Continuously  Stirred  Flow  Cell.  DP  Dissertation,  University  of  Florida, 
December  1994. 

Heyse,  Edward.  Faculty,  Department  of  Engineering  and  Environmental  Management, 
Air  Force  Institute  of  Technology,  Wright-Patterson  AFB,  OH,  Telephone  and 
Personal  Interview,  Summer  1995. 

Hildebrand,  F.  B.  Introduction  to  Numerical  Analysis  (Second  Edition).  New  York: 
Dover  Publications,  Inc.,  1987. 

Isaacson,  Eugene  and  Herbert  Bishop  Keller.  Analysis  of  Numerical  Methods.  New 
York:  Dover  Publications,  Inc.,  1994. 

Linz,  Peter.  Analytical  and  Numerical  Methods  for  Volterra  Equations.  Philadelphia: 
Siam,  1985. 

Masters,  Gilbert  M.  Introduction  to  Environmental  Engineering  and  Science.  Englewood 
Cliffs,  NJ:  Prentice  Hall,  1991. 

Microsoft  FORTRAN  PowerStation:  Professional  Development  System  Version  1.0  for 
MS-DOS  and  Windows  Operating  Systems.  Computer  Software.  Microsoft 
Corporation,  1993. 

Nkedi-Kizza,  et  al.  “Ion  Exchange  and  Diffusive  Mass  Transfer  During  Miscible 

Displacement  Through  an  Aggregated  Oxisol,”  Soil  Science  Society  of  America 
Journal.  46:  471-476(1982). 

- .  “On  the  Equivalence  of  Two  Conceptual  Models  for  Describing  Ion  Exchange 

During  Transport  Through  an  Aggregated  Oxisol,”  Water  Resources  Research. 
20181:  1123-1130(1984). 

Olsen,  Roger  L.  and  Michael  C.  Kavanaugh.  “Can  Groundwater  Restoration  Be 
Achieved.”  Water  Environment  and  Technology.  5:  42-47(1993). 


BIB.2 


Rao,  P.  S.  C.  et  al.  “Experimental  and  Mathematical  Description  of  Nonadsorbed  Solute 
Transfer  by  Diffusion  in  Spherical  Aggregates,”  Soil  Science  Society  of  America 
Journal.  44:  684-688(1980). 

Travis,  Curtis  C.  and  Carolyn  B.  Doty.  “Can  Contaminated  Aquifers  at  Superfund  Sites 
be  Remediated?,”  Environmental  Science  Technology.  24:  1464-1466  (1990). 

Treybal,  Robert  E.  Mass  Transfer  Operations  (Third  Edition),  New  York:  McGraw-Hill 
Book  Company,  1980. 

Tricomi,  F.  G.  Integral  Equations.  New  York:  Dover  Publications,  Inc.,  1985. 

Valocchi,  A.  J.  “Effect  of  Radial  Flow  on  Deviations  from  Local  Equilibrium  During 
Sorbing  Solute  Transport  Through  Homogeneous  Soils,”  Water  Resources 
Research.  22(12):  1693-1701  (1986) 

Vest,  Gary  D.,  Deputy  Assistant  Secretary  of  the  Air  Force  for  Environmental  Safety  and 
Occupational  Health.  Address  to  Air  Force  Institute  of  Technology  students  and 
faculty.  Air  Force  Institute  of  Technology,  Wright-Patterson  AFB,  OH,  23  July 
1992. 

Villermaux,  Jacques.  “Deformation  of  Chromatographic  Peaks  Under  the  Influence  of 
Mass  Transfer  Phenomena,”  Journal  of  Chromatographic  Science:  12:  822-83 1 
(1974). 

Weber,  Walter  J.  et  al.  “A  Distributed  Reactivity  Model  for  Sorption  by  Soils  and 

Sediments.  1.  Conceptual  Basis  and  Equilibrium  Assessments.”  Environmental 
Science  and  Technology.  26:  1955-1962(1992). 

- .  “Review  Paper:  Sorption  Phenomena  in  Subsurface  Systems:  Concepts,  Models 

and  Effects  on  Contaminant  Fate  and  Transport,”  Water  Res..  25(5):  499-528 
(1991). 

Wolfram,  Stephen.  Mathematica:  A  System  for  Doing  Mathematics  by  Computer 

(Second  Edition).  Reading,  MA:  Addison-Wesley  Publishing  Company,  1991. 

Wu,  Shian-chee  and  Philip  M.  Gschwend.  “Sorption  Kinetics  of  Hydrophobic  Organic 
Compounds  to  Natural  Sediments  and  Soils,”  Environmental  Science  and 
Technology.  20:  717-725  (1986). 

Yosida,  Kosaku.  Lectures  on  Differential  and  Integral  Equations.  New  York:  Dover 
Publications,  Inc.,  1991. 


BIB.3 


Young,  David  M.  and  Robert  Todd  Gregory.  A  Survey  of  Numerieal  Mathematics.  New 
York:  Dover  Publications,  Inc.,  1988. 


BIB.4 


Vita 


Capt  Jon  E.  Hodge  was  born  on  9  August  1963  in  Kansas  City,  Missouri.  He 
graduated  from  the  last  graduating  class  of  Olathe  High  School  in  1981,  and  entered  the 
Air  Force  Academy.  He  there  recieved  a  Bachelor  of  Science  degree  in  Physics  with 
two  tracks:  Optics  and  Aeronautical  Engineering,  Propulsion.  Upon  graduation  in  May 
1985,  he  entered  Undergraduate  Pilot  Training  at  Vance  AFB,  Oklahoma.  Graduating  a 
year  later,  he  underwent  six  months  of  weapons  system  specific  training  in  the  B-52G  at 
Castle  AFB,  California.  His  first  assignment  was  Blytheville  AFB,  Arkansas,  soon  to 
be  Eaker  AFB,  and  five  years  later  to  be  closed.  While  stationed  there,  he  flew  six 
combat  missions  from  the  island  of  Diego  Garcia  during  Desert  Storm.  He  then  began 
to  pursue  graduate  education  at  the  Air  Force  Institute  of  Technology  (AFIT).  After  two 
years  more  flying  in  the  B-52H  at  Minot  AFB,  North  Dakota,  including  a  demanding 
flight  scheduler’s  position,  he  entered  AFIT’s  School  of  Engineering  in  May  1994. 

Permanent  Address:  19615  K-68Hwy 

Paola,  KS  66071 


V.l 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  Instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collertion  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202*4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20503. 

1.  AGENCY  USE  OHIY  (Leave  blarik)  I  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

December  1995 _ _ Master's  Thesis _ _ 

4.  TITLE  AND  SUBTITLE  15.  FUNDING  NUMBERS 

A  POINT  MODEL  OF  AQUIFER  CLEANUP  WITH  A  DISTRIBUTION  OF 
FIRST-ORDER  RATE  PARAMETERS 

6.  AUTBOR(S) 

Jon  E.  Hodge,  Capt,  USAF 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES)  8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

Air  Force  Institute  of  Technology,  AFIT/ENP/GAP/95D-08 

WPAFB  OH  45433-6583 

9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES)  10.  SPONSORING /MONITORING 

AGENCY  REPORT  NUMBER 


11.  SUPPLEMENTARY  NOTES 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT  12b.  DISTRIBUTION  CODE 

Approved  for  public  release;  distribution  unlimited 

13,  ABSTRACT  (Maximum  200  words) 

Many  try  modeling  groundwater  contaminant  transport  to  predict  it.  Is  this  possible  with  rate-limited  processes,  and 
under  what  conditions?  On  occasion,  cleanups  go  slower  than  predicted  (tailing)  and  hazardous  concentrations 
reappear  after  cleanup  is  thought  complete  (rebound).  Rate-limited  transport  is  blamed  by  many.  When  immobile 
water  is  present,  diffusion  from  varied  sizes  and  shapes  of  immobile  regions  can  cause  varied  rate  limitations  (due  to 
varied  diffusion  path  lengths).  Although  known,  most  modelers  represent  these  varied  rate-limiting  processes  with  a 
single  “representative”  rate-parameter.  This  can  yield  poor  predictions  for  long-term  experiments,  and  the  parameter 
is  generally  time  and  pump-rate  dependent.  This  model  employs  a  distribution  of  first-order  rate  parameters  to 
investigate  the  effects  of  using  a  single  rate-parameter.  Spatial  effects  are  ignored  by  using  volume-averaged 
concentrations  (a  point,  well-mixed  model)  and  dilutive  pumping  and  rate-limited  transport  are  modeled  to  isolate  rate- 
limited  transport  for  study.  A  three-parameter  Gamma  distribution  defines  the  rate  parameter  continuum.  A  clean 
flow  approximation  is  used  extensively,  and  pulsed  pumping  is  examined  briefly.  An  effective  time  and  pump-rate 
dependence  is  seen  in  the  average  rate.  Long-term  soil  and  contaminant  transport  characteristics  along  with  uptake 
history  or  good  experimentally-derived  initial  contaminant  presence  are  concluded  as  necessary  for  accurate 
predictions. 

14.  SUBJECT  TERMS  15.  NUMBER  OF  PAGES 

Distributed  Model,  Gamma  Model,  Rate-limited  Sorption,  ,Tailing,  Rebound,  123 

Groundwater  Contaminant  Transport,  Aquifer  Decontamination,  Groundwater  Cleanup  16.  PRICE  CODE 

17.  SECURITY  CLASSIFICATION  I  18.  SECURITY  CLASSIFICATION  I  19.  SECURITY  CLASSIFICATION  20.1]mITATION 
OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 

Unclassified  Unclassified  Unclassified  UL 


NSN  7540-01-280-5500 


Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std.  Z39-18 
298-102 


GENERAL  INSTRUCTIONS  FOR  COMPLETING  SF  298 


The  Report  Documentation  Page  (RDP)  is  used  in  announcing  and  cataloging  reports.  It  is  important 
that  this  information  be  consistent  with  the  rest  of  the  report,  particularly  the  cover  and  title  page. 
Instructions  for  filling  in  each  block  of  the  form  follow  It  is  important  to  stay  within  the  iines  to  meet 
opticai  scanning  requirements. 


Block  1.  Agency  Use  OnW  (Leave  biank). 

Block  2.  Report  Date.  Full  publication  date 
including  day,  month,  and  year,  if  available  (e.g.  1 
Jan  88).  Must  cite  at  least  the  year. 

Blocks.  Type  of  Report  and  Dates  Covered. 
State  whether  report  is  interim,  final,  etc.  If 
applicable,  enter  inclusive  report  dates  (e.g.  1 0 
jun87-30Jun88). 

Block  4.  Title  and  Subtitle.  A  title  is  taken  from 
the  part  of  the  report  that  provides  the  most 
meaningful  and  complete  information.  When  a 
report  is  prepared  in  more  than  one  volume, 
repeat  the  primary  title,  add  volume  number,  and 
include  subtitle  for  the  specific  volume.  On 
classified  documents  enterthe  title  classification 
in  parentheses. 

Blocks.  Funding  Numbers.  To  include  contract 
and  grant  numbers;  may  include  program 
element  number(s),  project  number(s),  task 
numberfs),  and  work  unit  number(s).  Use  the 
following  labels; 

C  -  Contract  PR  -  Project 

G  -  Grant  TA  -  Task 

PE  -  Program  WU  -  Work  Unit 

Element  Accession  No. 

Blocks.  Author(s).  Namefs) of person(s) 
responsible  for  writing  the  report,  performing 
the  research,  or  credited  with  the  content  of  the 
report.  If  editor  or  compiler,  this  should  follow 
the  name(s). 

Block?.  Performing  Organization  Name(s)  and 
Address(es).  Self-explanatory. 

Block  8.  Performing  Organization  Report 
Number.  Enter  the  unique  alphanumeric  report 
number(s)  assigned  by  the  organization 
performing  the  report. 

Blocks.  Sponsoring/Monitoring  Agency  Name(s) 
and  Address(es).  Self-explanatory. 

Block  10.  Sponsoring/Monitoring  Agency 
Report  Number,  (if  known) 

Block  11.  Supplementary  Notes.  Enter 
information  not  included  elsewhere  such  as: 
Prepared  in  cooperation  with...;  Trans,  of...;  To  be 
published  in....  When  a  report  is  revised,  include 
a  statement  whether  the  new  report  supersedes 
or  supplements  the  older  report. 


Block  12a,  Distribution/ Availability  Statement. 
Denotes  public  availability  or  limitations.  Cite  any 
availability  to  the  public.  Enter  additional 
limitations  or  special  markings  in  all  capitals  (e.g. 
NOFORN,  REL,  ITAR). 

DOD  -  See  DoDD  5230.24,  "Distribution 
Statements  on  Technical 
Documents." 

DOE  -  See  authorities. 

NASA  -  See  Handbook  NHB  2200.2. 

NTIS  -  Leave  blank. 

Block  12b.  Distribution  Code. 

DOD  -  Leave  blank. 

DOE  -  Enter  DOE  distribution  categories 
from  the  Standard  Distribution  for 
Unclassified  Scientific  and  Technical 
Reports. 

NASA  -  Leave  blank. 

NTIS  -  Leave  blank. 

Block  13.  Abstract.  Include  a  brief  ('Max/mum 
200  words)  factual  summary  of  the  most 
significant  information  contained  in  the  report. 

Block  14.  Subject  Terms.  Keywords  or  phrases 
identifying  major  subjects  in  the  report. 

Block  1 5.  Number  of  Pages.  Enter  the  total 
number  of  pages. 

Block  16.  Price  Code.  Enter  appropriate  price 
code  (NTIS  only). 

Blocks  17.  - 19.  Security  Classifications.  Self- 
explanatory.  Enter  U.S.  Security  Classification  in 
accordance  with  U.S.  Security  Regulations  (i.e., 
UNCLASSIFIED).  If  form  contains  classified 
information,  stamp  classification  on  the  top  and 
bottom  of  the  page. 

Block  20.  Limitation  of  Abstract.  This  block  must 
be  completed  to  assign  a  limitation  to  the 
abstract.  Enter  either  UL  (unlimited)  or  SAR  (same 
as  report).  An  entry  in  this  block  is  necessary  if 
the  abstract  is  to  be  limited.  If  blank,  the  abstract 
is  assumed  to  be  unlimited. 


★  U.S.GPO:  1 993-0-336-043 


Standar(d  Form  298  Back  (Rev.  2-89) 


